{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analysis to calculate the diffusivities, Nernst-Einstein conductivites, radial distribution function, coordination number and ion pair lifetimes of molten NaCl\n",
    "\n",
    "Assumes the path to PyLAT/src has been added to PYTHONPATH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import Classes\n",
    "\n",
    "import and initialize required classes for calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MSD import MSD\n",
    "from calcDiffusivity import calcdiffusivity\n",
    "from calcCOM import calcCOM\n",
    "from getTimeData import gettimedata\n",
    "from getMolData import getmoldata\n",
    "from COMradial import COMradialdistribution\n",
    "from getAtomCharges import getatomcharges\n",
    "from calcNEconductivity import calcNEconductivity\n",
    "from getCoordinationNumber import getcoordinationnumber\n",
    "from ionpair import ionpair\n",
    "\n",
    "\n",
    "c = calcCOM()\n",
    "m = MSD()\n",
    "cd = calcdiffusivity()\n",
    "gt = gettimedata()\n",
    "gm = getmoldata()\n",
    "crd = COMradialdistribution()\n",
    "gc = getatomcharges()\n",
    "ne = calcNEconductivity()\n",
    "ip = ionpair()\n",
    "gc = getatomcharges()\n",
    "gcn = getcoordinationnumber()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# General Input\n",
    "input used by all calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = {}\n",
    "datfilename = 'mol.data'\n",
    "trjfilename=['mol.lammpstrj']\n",
    "logfilename='log.lammps'\n",
    "temp=1300\n",
    "nummoltype=[864,864]\n",
    "moltypel=['Na','Cl']\n",
    "moltype = []\n",
    "for i in range(0,len(moltypel)):\n",
    "    for j in range(0,nummoltype[i]):\n",
    "        moltype.append(int(i))\n",
    "verb=2\n",
    "ver=True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Input for Radial Distribution Function and Coordination Number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "RDF_maxr=None\n",
    "RDF_binsize=0.1\n",
    "RDF_Timesteps = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Input for Ion Pair Lifetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "IPL_skip=0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Input for Diffusivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "MSD_skip=0\n",
    "MSD_num_init=None\n",
    "D_Tolerance=0.075"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preliminary Calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "COM calculation 100.00% complete\n"
     ]
    }
   ],
   "source": [
    "output['Conductivity'] = {}\n",
    "output['Conductivity']['units'] = 'S/m'\n",
    "\n",
    "n = gc.findnumatoms(datfilename)\n",
    "(molcharges, atomcharges,n) = gc.getmolcharges(datfilename,n)\n",
    "molcharge = gc.molchargedict(molcharges, moltypel, moltype)\n",
    "tsjump = gt.getjump(trjfilename[0])\n",
    "dt = gt.getdt(logfilename)\n",
    "(V, Lx, Ly, Lz) = gcn.getvolume(trjfilename[0])\n",
    "(comx, comy, comz, Lx, Ly, Lz, Lx2, Ly2, Lz2) = c.calcCOM(trjfilename,datfilename, ver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation of Radial Distribution Function and Coordination Number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "output['RDF'] = {}\n",
    "output['RDF']['units'] = 'unitless, angstroms'\n",
    "output = crd.runradial(datfilename, comx, comy, comz, Lx, Ly, Lz, Lx2, Ly2, Lz2, output, nummoltype, moltypel, moltype,RDF_Timesteps, ver,RDF_maxr,RDF_binsize)\n",
    "output = gcn.calccoordinationnumber(output, nummoltype, moltypel, V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting RDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f31b9a2b2d0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEQCAYAAABfiGi4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VNXZwPHfM5N9Yw2BECCACwjKFkEUFNxFBK2oULW4FRdstXWp+lK11tpqXVq1VXEr7rsIClVUqCsoIPsua2QLCWTfZua8f9xJTMJMMpmZOxPC8/Uzzsy95955mCTzzFnuOWKMQSmllAJwRDsApZRSLYcmBaWUUrU0KSillKqlSUEppVQtTQpKKaVqaVJQSilVy/akICJOEflBRD70sS9eRN4UkU0iskhEsu2ORymllH+RqCncBKz1s+9qYL8x5gjgMeDBCMSjlFLKD1uTgohkAecCz/kpMh6Y4X38DnCaiIidMSmllPLP7prCP4DbAY+f/V2BHQDGGBdQCHSwOSallFJ+xNh1YhEZC+w1xiwRkVEhnmsKMAUgOTl5SJ8+fcIQoVJKHT6WLFmyzxiT3lQ525ICcBIwTkTGAAlAmoi8Yoy5rE6Zn4BuQK6IxABtgPyGJzLGTAemA+Tk5JjFixfbGLZSSrU+IrItkHK2NR8ZY+40xmQZY7KBicDnDRICwCxgsvfxBG8ZnaFPKaWixM6agk8ich+w2BgzC3geeFlENgEFWMlDKaVUlEQkKRhjFgALvI/vrrO9ArgoEjEopZRqWsRrCkopZYfq6mpyc3OpqKiIdihRlZCQQFZWFrGxsUEdr0lBKdUq5ObmkpqaSnZ2Nofr5U7GGPLz88nNzaVnz55BnUPnPlJKtQoVFRV06NDhsE0IACJChw4dQqotaVJQSrUah3NCqBHqe6BJoSXYtRy2L4p2FEoppUmhRZh3D3x4c7SjUEqFSES45ZZbap8//PDD3Hvvvc06x6hRo8jJyal9vnjxYkaNGhWmCJumSaElKNkL+7eBXren1CEtPj6e9957j3379oV0nr179zJ37twwRdU8mhRagtI8qC6FsoJoR6KUCkFMTAxTpkzhscceO2jf7NmzGTZsGIMGDeL0009nz549fs9z22238Ze//OWg7Vu3bmXkyJEMHjyYwYMH880334Q1ftAhqdHn8UCZ91tF4XZI1klilQrVn2avZs3OorCe85jMNO45r1+T5aZOncpxxx3H7bffXm/7iBEjWLhwISLCc889x0MPPcQjjzzi8xzDhw/n/fffZ/78+aSmptZu79SpE/PmzSMhIYGNGzcyadIkwj0XnCaFaCvfD8Y7s/iBHZA5KLrxKKVCkpaWxq9+9Ssef/xxEhMTa7fn5uZyySWXsGvXLqqqqpq8jmDatGncf//9PPjgz2uPVVdXc+ONN7Js2TKcTicbNmwIe/yaFKKtNO/nxwe2Ry8OpVqRQL7R2+nmm29m8ODBXHnllbXbfvOb3/D73/+ecePGsWDBgtoO6LPOOos9e/aQk5PDc8/9vB7ZqaeeyrRp01i4cGHttscee4yMjAyWL1+Ox+MhISEh7LFrUoi2ukmhcEf04lBKhU379u25+OKLef7557nqqqsAKCwspGvXrgDMmDGjtuzHH3/s9zzTpk3juuuuo1evXrXnyMrKwuFwMGPGDNxud9hj147maKtJCs44rSko1Yrccsst9UYh3XvvvVx00UUMGTKEjh07BnSOMWPGkJ7+87o4N9xwAzNmzGDAgAGsW7eO5OTksMcth9ryBa1ukZ1F02HubdB1CLiq4Pqvoh2RUoektWvX0rdv32iH0SL4ei9EZIkxJsfPIbW0phBtpXmAQJcB1ugjpZSKIk0K0VaaB0kdoF02VBRaN6WUihJNCtFWmgcpnaBtd+v5Ae1sVkpFj21JQUQSROQ7EVkuIqtF5E8+ylwhInkissx7u8aueFqs0n2Q3BHa1CQFbUJSSkWPnUNSK4FTjTElIhILfCUic40xCxuUe9MYc6ONcbRspXmQOfDnK5nLdaoLpVT02JYUjDWsqcT7NNZ7O7SGOkVC6T5ITof4NOt5ZXF041FKHdZs7VMQEaeILAP2AvOMMb4WDbhQRFaIyDsi0s3OeFocVyVUFlrNR/He+U0qwjtfi1IqcsIxdXZ1dTV33HEHRx55JIMHD2b48OG1M6ZmZ2eHPANrU2xNCsYYtzFmIJAFDBWR/g2KzAayjTHHAfOAGQ3PASAiU0RksYgszsvL81Xk0FTq/eEmp4MzFmKToFKTglKHqnBMnf3HP/6RXbt2sWrVKpYuXcrMmTMpLo5cC0JERh8ZYw4A84GzG2zPN8ZUep8+Bwzxc/x0Y0yOMSan7tV9h7yaq5mTvf+m+FRNCkodwkKdOrusrIxnn32WJ554gvj4eAAyMjK4+OKLbY+9hm19CiKSDlQbYw6ISCJwBvBggzJdjDG7vE/HAWvtiqdFqjhg3Se2s+7j07RPQalwmHsH7F4Z3nN2PhbO+VuTxUKZOnvTpk10796dtLS0sIbeHHaOPuoCzBARJ1aN5C1jzIcich+w2BgzC/itiIwDXEABcIWN8bQ8Lm8lKcY702F8qvYpKHWIC9fU2dFi5+ijFcBBiwMYY+6u8/hO4E67YmjxXBXWfU1SSNCaglJhEcA3ejsFO3X2448/zvbt2ykqKopabUGvaI6m2pqC1XaofQpKtQ51p86u0djU2cuWLeO5554jKSmJq6++mptuuomqqioA8vLyePvttyMWuyaFaGpYU4hvozUFpVqJYKfOvv/++0lPT+eYY46hf//+jB07NqK1Bl1kJ5q0T0GpVqWkpKT2cUZGBmVlZbXPx48fz/jx45s8R1xcHA899BAPPfTQQfu2bt0aljgbozWFaKqtKXibjxLSoKoYPJ7oxaSUOqxpUoimg5qPvFc1V2kTklIqOjQpRJOrEsQJTm8rns5/pJSKMk0K0eSq+LmWADr/kVIq6jQpRJOr8uf+BLD6FEBrCkqpqNGkEE0H1RRqkoLWFJRS0aFJIZqqK+rXFDQpKHXI2717NxMnTqR3794MGTKEMWPGsGHDBvr3bzhJ9M8efvhh+vTpw8CBAzn++ON56aWXABg1ahSLFy+OVOiAXqcQXdqnoFSrYozhggsuYPLkybzxxhsALF++3OeMqDWefvpp5s2bx3fffUdaWhpFRUW8//77kQr5IJoUokn7FJRqVebPn09sbCzXXXdd7bYBAwY0etHZAw88wIIFC2qvWk5LS2Py5Ml2h+qXJoVoalhTiE0GRJuPlArRg989yLqCdWE9Z5/2ffjD0D80WmbVqlUMGeJzWRifioqKKC4uplevXqGGFzbapxBNDWsKDoeuqaCUiiqtKUSTqwIS2tTfpvMfKRWypr7R26Vfv3688847jZa58sor+eGHH8jMzGTOnDmkpKSwefPmFlNb0JpCNDWsKYB3TQVNCkodik499VQqKyuZPn167bYVK1awY8eO2ucvvvgiy5YtY86cOQDceeedTJ06laIi6+++pKSkdvRRNGhSiKaGfQqgayoodQgTEd5//30+/fRTevfuTb9+/bjzzjvp3Lmz32Ouv/56Ro8ezfHHH0///v0ZOXIkDkf0PprtXKM5AfgCiPe+zjvGmHsalIkHXgKGAPnAJcaYrXbF1OL4qinEp0FpXnTiUUqFLDMzk7feeuug7atWrfJZXkS4/fbbD1rTGWDBggXhDq9JdqajSuBUY8wAYCBwtoic0KDM1cB+Y8wRwGPAgzbG0/L4rSloR7NSKjpsSwrGUrPiRKz3ZhoUGw/UrEv3DnCaiIhdMbU4PmsKmhSUUtFja8OViDhFZBmwF5hnjFnUoEhXYAeAMcYFFAId7IypRfFXU6gq8V1eKdUoYxp+7zz8hPoe2JoUjDFuY8xAIAsYKiL+J/9ohIhMEZHFIrI4L6+VtLe7XWDcvpNCdRl43NGJS6lDVEJCAvn5+Yd1YjDGkJ+fT0JCQtOF/YjIdQrGmAMiMh84G6jb2/IT0A3IFZEYoA1Wh3PD46cD0wFycnJax0+84VKcNeJSrPvKYkhsG9mYlDqEZWVlkZubS6v54hikhIQEsrKygj7eztFH6UC1NyEkAmdwcEfyLGAy8C0wAfjcHC5p3lVp3R9UU/AmhaoSTQpKNUNsbCw9e/aMdhiHPDtrCl2AGSLixGqmessY86GI3AcsNsbMAp4HXhaRTUABMNHGeFqWJmsK2q+glIo825KCMWYFMMjH9rvrPK4ALrIrhhatNik0rCl4Z0rVzmalVBToFc3RUtt81HBIap0+BaWUijBNCtHir6YQp0lBKRU9mhSipamagjYfKaWiQJNCtDTVp6AdzUqpKNCkEC3+ago1zUdV2nyklIo8TQrR4q+mEBMPjhitKSilokKTQrT4qymIWLUF7WhWSkWBJoVo8VdTAJ0UTykVNZoUosXfNBeg02crpaJGk0K0+JvmAqzmI60pKKWiQJNCtDRaU0jRjmalVFRoUogWVwWIE5w+pp/SjmalVJRoUogWX6uu1YhP0+YjpVRUaFKIFl/rM9fQ5iOlVJRoUoiWxmoKcSnWFc2HyXpDSqmWQ5NCtDRVUzAeqC6PbExKqcOeJoVoabRPIdW6185mpVSE2ZYURKSbiMwXkTUislpEbvJRZpSIFIrIMu/tbl/napUaqynEeZOCdjYrpSLMzjWaXcAtxpilIpIKLBGRecaYNQ3KfWmMGWtjHC1TozUFXWhHKRUdttUUjDG7jDFLvY+LgbVAV7te75DTaE1BF9pRSkVHRPoURCQbGAQs8rF7uIgsF5G5ItIvEvG0CK6KRjqaa/oUNCkopSLLzuYjAEQkBXgXuNkYU9Rg91KghzGmRETGADOBI32cYwowBaB79+42Rxwh7ipwxvnepx3NSqkosbWmICKxWAnhVWPMew33G2OKjDEl3sdzgFgR6eij3HRjTI4xJic9Pd3OkCOnsaSgq68ppaLEztFHAjwPrDXGPOqnTGdvOURkqDeefLtialHcVQE0H2lSUEpFlp3NRycBlwMrRWSZd9tdQHcAY8zTwATgehFxAeXARGMOk8t43dXgjPW9Ly7ZmiyvomFrm1JK2cu2pGCM+QqQJso8CTxpVwwtWmPNRyLehXY0KSilIkuvaI4Wd7X/pACQkKY1BaVUxGlSiBZXpf/mI4D4NlpTUEpFnCaFaDCm8eYj0JqCUioqNClEg8cNGHD6GX0E1kI7lYURC0kppUCTQnS4q6z7xpqPtKaglIoCTQrRUJsUGmk+ik/TPgWlVMRpUogGd7V1H0hN4TC5bEMp1TJoUogGd6V131RNwbihqjQyMSmlFJoUoiOQ5qOENOtem5CUUhGkSSEaapqPYpqoKYB2NiulIkqTQjQEVFNoY91rTUEpFUGaFKIh0NFHoDUFpVREaVKIhkBHH4FewKaUiihNCtHgCnD0EWhNQSkVUZoUoqG2ptDINBc6+kgpFQWaFKIhkGku4lJAHFpTUEpFlCaFaAiko1lEp7pQSkWcnWs0dxOR+SKyRkRWi8hNPsqIiDwuIptEZIWIDLYrnhaltvmokaQAOimeUiri7Fyj2QXcYoxZKiKpwBIRmWeMWVOnzDnAkd7bMOAp733rFkjzEehCO0qpiAsoKYhIJ+AkIBMoB1YBi40xHn/HGGN2Abu8j4tFZC3QFaibFMYDLxljDLBQRNqKSBfvsa1XIHMfgdYUlFIR12hSEJHRwB1Ae+AHYC+QAJwP9BaRd4BHjDGNfnKJSDYwCFjUYFdXYEed57neba08KdRMc9HI6COw+hSKcu2PRymlvJqqKYwBfm2M2d5wh4jEAGOBM4B3/Z1ARFK8+29uKnk0co4pwBSA7t27B3OKliXQ5qOENNirNQWlVOQ02tFsjLkNyBWRi33scxljZhpjGksIsVgJ4VVjzHs+ivwEdKvzPMu7reFrTTfG5BhjctLT0xsL+dAQyOgj0NFHSqmIa3L0kbff4PbmnlhEBHgeWGuMedRPsVnAr7yjkE4AClt9fwL83HzkaKqm0EYX2lFKRVSgo48+FZFbgTeB2lVfjDEFjRxzEnA5sFJElnm33QV09x77NDAHq4lqE1AGXNms6A9V7ipwxICjiZyc2M5aaKeiEBLbRiY2pdRhLdCkcIn3fmqdbQbo5e8AY8xXgDR2Uu+oo6mNlWmV3FWNT3FRI6m9dV9eoElBKRURASUFY0xPuwM5rLiqmu5kBkjqYN2X7bfGfymllM0abb8QkRFN7E8Tkf7hDekw4K5qupMZINGbCcry7Y1HKaW8mqopXCgiDwH/BZYAeVjXKRwBjAZ6ALfYGmFr5K4OLCnUbT5SSqkIaDQpGGN+JyLtgQuBi4DOWFc0rwWeNsZ8bX+IrZA7wOajxHbWfZkmBaVUZDTZp2CMKRCRNGAFsLJmM9BHREqNMcv8H618CrT5KKGtNX221hSUUhES6CypQ4DrgC5Y8x9dC5wNPCsizb6G4bDnroaYAJKCw2ElBq0pKKUiJNAhqVnAYGNMCYCI3AN8BJyM1dfwkD3htVLuysBqCmD1K2hNQSkVIYHWFDoBlXWeVwMZxpjyBttVIAJtPgJrBJLWFJRSERJoTeFVYJGIfOB9fh7wmogkU38qbBUId3VgHc1g1RSKDpoOSimlbBHoxWt/FpG5WFNXAFxnjFnsfXypLZG1Zu4qiEsOrGxie9i9yt54lFLKK+CV17xJYHGTBVXTmtN8pH0KSqkIsm2NZtWIQC9eA+taheoyqC63NyallEKTQnS4mjn6CELubN5csIeK6qqQzqGUav00KURDs2oKoU11UVFdxVUz/8q4WWcx4uUJ/LBzS1DnUUodHjQpREOg01xAnZlSm58UCivKOP3VK/i+8DXa0JcK2cmv/juJL7fogDGllG+aFKKhuR3N0OyaQpXLxTmvX8kBVnFmxg18fcXr/HvUDABu/PxGth/Ia9b5lFKHB00K0RDoNBdQZ/rs5iWFJxfOotixhjMzruWRs68H4OSe/bhryIO4HQX8cuZNeDyeZp1TKdX62ZYUROQFEdkrIj4H2YvIKBEpFJFl3tvddsXS4kSgpvDWxtcQVzseOP3X9bZPGnAKp3a6mkJZyR/mPduscyqlWj87awr/wZo0rzFfGmMGem/32RhLy+HxgKcZHc0x8RCXCqX7An6Jjzf+QKljPSekn0dC7MGv8+jZN5Di6cvcndNZtmtrwOdVSrV+tiUFY8wXgF511ZCn2roPtKMZIDUDincHXPzfi1/BeGK4a8Rkn/tjnE6eOOMBwMPtnz0YeBxKqVYv2n0Kw0VkuYjMFZF+UY4lMtzeawUCrSkApHSGkj0BF99atpR2jv5kt+/kt0xO1hEckzKGna6vmb32+8BjUUq1atFMCkuBHsaYAcATwEx/BUVkiogsFpHFeXmH+KgZd01NoRlJoRk1he9zN+GJ2cfAjsc3WfbRM29BPIk8sPDv2umslAKimBSMMUU16zMYY+YAsSLS0U/Z6caYHGNMTnp6ekTjDLtQagrGNFn0ndXzAbigz6gmy2a1ac8pnX5JiWMtT333UeDxKKVaraglBRHpLCLifTzUG0t+tOKJmGCSQmqGNf9RZXGTRRfv/Q7caYzq2T+gUz94xnU4XB15bvUTOg2GUsrWIamvA98CR4tIrohcLSLXich13iITgFUishx4HJhoTABfhQ91riBrCtBkv4LL7WZv9Sq6xB6LwxHYjzY5Pp5JR16HK2YXf/vyjcBjUkq1SgFPnd1cxphJTex/EnjSrtdvsWprCs0cfQRWv0LHI/0Wm795JThLOL7z0GaFdOtJE3hz4wt8sOUVprkvJcbpbNbxSqnWI9qjjw4/wfYpQJM1hQVblwJwRu/mJYUYp5MLel2OK2YXj3/7QdMHKKVaLU0KkVYz+ijQaS6gfk2hEav2rcZ44hnRo2+zw7p9xMU4XB14bcOLOhJJqcOYJoVIC6amkNAWnPFQ0nhS2Fm+iSTTPajmn4TYOM7MmkSlcysvLp3X7OOVUq2DJoVICyYpiHivVfDffFRRXUW57CAryX+fQ1PuHvUrcKfx7Mrngj7HQYyx4va4w3dOpZRtNClEWjAdzeC9VsF/TeHLrWsQRzXHdQr+wvDU+ERGpk+g1LGOd1Z+HfR5AGuupplT4a9Z8MhR8PBR8NGtUNr6Rx0rdSjTpBBpwdQUoMmawv+2/QDAKdmDgo0MgD+Nvho8CTyz7KXgT7L1K3gyB1a8Cf1/AWf9FXqdAktetLav/29IMSql7GPbkFTlR7BJIaUzbPnC7+5V+9ZgPHGc1L35ncx1paekcUTiaDZWfML6vJ0cnZ7ZvBOsmwNvXwHtsuGqjyH96J/37VkDM6+DN34J4x6HQZeFFKtSKvy0phBprkrrPia+ecelZkBFIVSX+9y9s+xHkkw34mJCz/O/H3YFIm7+9tV/mndgwRZ479eQcQxc9d/6CQGs7VfMgZ4nwwdTYeU7IceqlAovTQqRVvOhHpPYvOPadLPuC3MP2uXxeCgnl04J2aHF5jWy5zGkevqxZP+cwKe+cLvg/WtBnHDxyz8vDtRQfApMegO6nwgzr4etIfZdKKXCSpNCpLkqrPvYhOYd1y7but+/9aBda/bmgrOcI9oFP/KooYuPnohxFvL4t34nr61v6X9gxyI492Fo263xsrEJMPFVaNsD3p4MRTtDjlcpFR6aFCKttqbQzKTQtod17yMpfLndWvF0UEafEAKr74ahYxFXe9798c2mC1cWw4K/QY+T4NiLAnuBpPZWYqgqg7ev/PmiPqVUVGlSiDRXJSBBdDRnWInER1JYsWcdAKdkHxt6fF5xMTGckD6WMscG5m1c1njhb56A0jw448/WNRWBSj/a6nDesRA+vTekeMFqRttZVMCXW9awYPMqth84xNfeUCoKdPRRpLnKITaxeR+eAA6HVVs4sO2gXZsLN4E7rdGV1oJx14jJjP3gNf75/X8448h/+C5UvNtKCv0ugKwhzX+RYyfA9oXw7ZPQbRgcM67ZpyirruT++a/w0Y5X8cTUTwROVwZdEvowPPMEbhr+C9okJDU/RqUOI5oUIq26ovkjj2q06+GzprCvajupkhVaXD5kt+9EVuyJbK38kp1FBWSm+eg8XvBXq+nntLuDf6Gz/gI/LYFZN0LXwdAmsH+Lx+Nh2mcv8uH2lzAxBcRJN3LaTaZLSicc4mBb4U9sOLCS3MrveXv7/3h7yxOc2PEiHj3rNyTHB/kzUKqV06QQaa6K5o88qtEu2/pWbUxtTaPK5aJSdpKddGb4Yqzj+sG/Ytr3X/DXL1/hiXN/W39n3npY+hIMnQLtewX/IjHxMOF5eHokvHctTJ4Fjsbnb9pTUsik924hzywiXrK54uhbuWHouT7XkXC53cz44VNeXDWDb/bP4KRXPuaPJ/yJC/udGHzMSrVSmhQizVXR/JFHNdplQ2URlO+vHfK5+KdNiKOao8M48qiu8ccM488Le/Ll7g9wuafWn2zvq8esfo6Tbwv9hdr3gnMegg9usGofp07zW/TTTcu59X+34HLuZVjby3jmvFsbnQQwxunk6pyzuDrnLB7/5gOeW/cw93x/HW+vGc+z4+4iNT6IJF20CzbMhf3boGyfNa1H6T5wV0KnfpA5CLoOsWo+TSQ4pVoSTQqRVl3e/JFHNeqOQPImhUW5awEY2OVoPweF7rzsi3hnx0M8v+Rjrh06xtpYmAsr34bjfw3JPpfWbr6Bv4Tt38IXf7c+VPuce1CRN1d8yZ+X/A6RWG459hGuHHJGs17ityeOZ0L/kVw9+25Wl83k5FcWcc/w+zj/mBOaPriy2Jq6Y9V7sO0bwIAj1vr3J3W07h3tYPN8WOFdxS4lA445H/pfCFnHW31DSrVgtiUFEXkBGAvsNcYctGCwd33mfwJjgDLgCmPMUrviaTFclcEnhZprFQ5ss76BAuvyfwTgxO7HhCE4335/0gTeefXfvLL29Z+TwsKnrGas4TeE74VEYMzDsGeV1Yw0ZX69leZmrlnI/Ut+j9OTxsvnvsBxnbODepnMtPbMvfRJ/rVwNs+seZBp313L22svYPrYO333NVSWWJ3pC5+CykJI7wuj7oR+50PHo3wPGijaBdu+htXvw5L/wHfPQFoWDLgETpgKyR0CjrfK5WL2+u+Yt/lbdpfupsJVhsFQVL2PalNBWkw6HRM6k5WaSa923TkmvQeDuvSiXVJKUO+POrzZWVP4D9Zym/5mVjsHONJ7GwY85b1v3VwV1uijYLQ7+FqF7cXbwJ1CVhs/VxCHQWp8IsemncWKkvf4PncTx3foYH3Q9b8Q2nYP74vFJlhXRE8/Bd64FH79GcSnsnD7ev648HeISWLGmOATQl1TTziP8X1P4urZ01hR8i4jXvmW+068n/P6Hm8VMAbWzoK5f4DiXdB3HJx0E2TlNH3ytC7WyKpjJ0BFEayfC6vehS8fhUXPwIBJMGQydPY9jHjF7q28seIzFu7+hjzXKnCWWTs8CYhJAANx0o5YSSC/egt73UtYW+6GvcB6MEaI83QlM+FohmQM4qwjhnFCt6MCW7u7eI9VEyzNs261zWN51txdscnW73BckvU4LdP6PWjXA9K6Nn8GYNWiiDHGvpOLZAMf+qkpPAMsMMa87n2+HhhljNnV2DlzcnLM4sWLbYg2Qp49DRLS4PL3gzv+od5Ws8q4xwEY+sIvAMN3VwV5vgAt27WVyz4ez1EJZ/NeZjp8/me47mvofNCPNjy2fAEvjYc+Y9lx5mOMfe9iPFLKv0e/yMie4a8VPf7NBzy77iGMo4RBqRN4+pTLSP50Gqz70PrgHvMIdA/Dd5a96+CrR2H1TKv/IXMwDLkC+l9Ivttw9+fP882ej3DFeKdJd6eRGTeAk7qeyKXHnU7vDp19ntbldrN+305W7tnChvztbCzYzI9FqylmMzi8V9G7U2jnPJIBHXP4Rd/RnJLdD4erwhr59dMS+Gkx5C6BYh9XmMelQFIH6/qa6nKoLrUuPHRX1i8nTmjT1eojat8bOvSGTn2h0zFWU5oIHo+Hjfm7WbV3Kxvzd7Ct8Cd2le4iv3I3Ja69eKi5kNF4b2C894LglCTiJYVEZypJMamkxaXRNr4d7RIcZqGLAAAa9klEQVTa0DmlA0d17M6xGT3qj5YzxuqPK8u3+uTKD0DFgZ/vXVXgjLGma3FXWiPqHDFW8otJqH8flwLxqdbfcbz3lpDmNxl6PB6KKsvZU1JIfllR7a2kupx4Zxwu46a4ogSnw0l8TBzxzjgSYmKJc8aSEBNHQmwcHRJT6Z/RnYTYZl7fVPdHI7LEGNPkN5poJoUPgb8ZY77yPv8M+IMxptFP/EM+KTw1wvpWNem14I5/8Vzr29o11upox74wnKz4HOZe+kQYg/TtzFeuZ1f1Ij7fV0h6l4Fw2bv2vuA3T+D5ZBoTug5mQ+x+/m/QP5k04BTbXm77gTyunv1Hdnu+5qhKF/cXFNJ35B+s5h5nmCvVZQVW/8SSGaw/sIm30toyKzmZCqeHRHdvhqSfzPl9RnNG7wGBfbv3o8rlYv6WlXy6+TtW5C1nT+Ua3DHWmhbtXMKI8lKGVpRzVFUV2cldSeqaY9WE2vey+kiS063+kjg/13e4q6HoJziw3ep0P7Ad9m+Bgs2UFPzIZk8ZG2Nj2RQXy4a4RLbHxLHPCS5H/c8d44kn1tOBlJhOxDms5lXx/of3/4hgjKHcXUylp4RqU4JbSjGOMkQO/hyL9wgdXUJHj5t27kraeVykejykeTykuq37urcUj8GBwe2Iwe2IxXjcGFONQfAA1SIUORwUORwUOq37YoeDUodQ6nBQ7IihyBFDkTgpcTgocwhlYqhwGDzNvCzJF4cRTowbyVO//FdQxweaFA6JjmYRmQJMAejePczNFZHmKg/+OgWAjH6w7FXweMgtPgDOErqn9ghffI2YdtKNTP3fV7wc5+L3J91k/wsOv5EHVr3Hxrh9TIo9wdaEANA9MZF5KcLnG/O4u0M6F3dOZ9CeYv7lqiY1zEmhKi6Nf1R14K3ETlSmlOI0cHppGZcVHWBg2zRIc0FaXPMvcmwgzlPJWY4DnBWbB1VbYOdqfhI33yQm8VlyR+YlpTI71frAN8YQU7iVtqUeuqXs45iOR9G/k5usNm5S4hJIjU8gLSGJBGcsRZXl5JVa33gLykvYlL+fZXtzySvfx4HKCgodbkznNCANAKfHQRdXHL1c1ZxaXk736hIyXW4yXS66uFykOeIhtRriSiAu2brFJh387dvjgvJCqNxvJdbyAjyuCkpFOOB0kO90sjsmhp0xseyIS2GHM558Zyy7nImUiocKhxu3I7yrAMZ4nMR5HCQYB8lGSDbQyWNIdbtJ87hJ87ho464kzVVJivGQ7DGkeDzEG0O1CDHGkOgxGIFqhGoRXGIlIes5FDsc7IiNoUt6aVhj9/nvsf0V/PsJqDtzWpZ320GMMdOB6WDVFOwPzUbVIfQpgDX9dFUJFG5n0U6rmn90KNcINMPJPfpwWrmH19LSOC/lCOwZBPuzv3/1Lm/G7+Pk8gTu3PI2rDzLaqO3w47vrVle92/h1JNvo9exV3LN3Hv4ofhNRr4yn8l9buSm4eND+tbu8Xj4aP0SZqx8jw0lX2KchThMR07ucDW3nfhLspPirRFdy9+A+Q/A/L9YtcoeI6yBBR2PsqYG8TbD1ON2WSvzFe2CvWtg13LYtcy697isZp0uA2DYtXTtMYKLup/ARYltqXK5+Gb7OhbmrmZ13nq2l/zIftc29hUvYVmJga3N/Ee6k4gxbUmPO4oeqb04Nr0PJ3brz5CuveoPG64qtWoVtTWMbVCyx2qSqiqx+mGKdoFp8AEuDkhsZ70vXQZCUjscie1JTWpPakpnuqV1gdRMq6nLz8+qtLKS3SX72VlUwJ6S/eSVF5JfdoADFUUYDE5x4hDBIQ4c4sApDpwOJzEOJx2T2tIpuR1dUjqQmdaejJQ2gU9X7/FAVbH1b6tZV8V4rKYtTP3HcSlWc5XH5W26ircGqUSgvyaaSWEWcKOIvIHVwVzYVH9Cq+CqCH70EUCGtyVuz2pW7N4HwODMo8IQWADWzuLm/N18ltWV2z57jJkT/27bS/13w1JmbPorCaYnD1z6Ao63L4f3plh/ID6Gqgatuhzm3QPfTbc6TCfPhuwRZAOfXj6dfy/6kGdWPcoLm+7mpXVP07/tCK7LmcBJPQJfzGhx7iaeWvwOS/M/wxWzG2MctHMex3k9z+PmE39R/0Pl+GusW/Ee6zqIDZ/Axo9heZ3mRnFYQ2EdMdZNxFprgzrfl+LTrCQw/EbIHmn1h8SnHhRbXEwMo3r1Z1Sv+i28+8tK+GrbWtbt20p++QEq3FVUuCqodFVS5akiwZlISlwyqXHJtIlPITMtndN7DyIjpU1gb0pcsrevIbRFoYKRHB9P7/jOfvtnbONwQEIb6xaMuOTwxuOHbX0KIvI6MAroCOwB7gFiAYwxT3uHpD4JnI01JPXKpvoToBX0KTzQFQZPhrMfCO74yhJr3ePRd/HLvDJWlLzPoknf2T9tgzEwfRRUFnNuhxFsq/6C5099m2Hdw19f2LhvFxd+MBHEzZtj36Bvpyzr29XL51vffM96wLqKOsSmFfastoa+7lkJQ6+F0/7o84OzorqK++a/zKe5H1Lm+BERQ4K7F92T+zIo4zhO7NafjJR2VLld7CstZF9ZIcv3bGDlvlXsKt9AdYy1Bkai+whGZp7Jb4Ze2Lx5qoyx5pjat966irxkr/UN2uOyaggel3XdSmoX65Z+FLTN1msiVD0toqPZDod8UvhTexhxc2hzBf1zIHQ5jtEVTgqqt7H86s/DF58/mxdYo4HG/oPV3U7nkjnj6RyTw6eXTw/ry5RWVjL61UmUyRbuH/ZU/YvKKgqtD/ENc+Hoc2HMQwHPk1RPWYHVPLP4eUhoCxc8A0cFNk3I6j07eOTbV1hR8C0VsgNxuPwXdieSKj3p224w1+dcSE7WEc2PVakwaVUdza2G22V9wwt27qMaGf1gzxoOJHUk1dklPLE15et/QnInGDCJfrEJHJc6lpWl7zN77fc/j+sPkcfj4aJ3bqXcuZGLu//h4KuME9rAxNfg2ydg/l/hyaEw+i4Ydl1go4Mqi+H7563pOSqLrWaaUXf6XyXOh34Z3Xjh/DsBa3bWzzat4Ifd6ymsLCbWEUOb+DTaJCRzbKdeDO9+dKPTbyjVEmlSiCRXzQI7ITb1ZPTDvX4OrhQPGYkDQ4+rKbuWw4+fw2n31M7b9PCZN3PWOx/z14WPcF7fN8LyMtfMepAdrgUcm3wBd4++zHchh8O6gOyY82HObfDJ/1lXGx87AQZMrH8xmDHWqm4/fg6b5sGmz62Ovt6nwZn3W532IUiKjee8vseHLSkq1RJoUoik6pqlOEOsKXQ+jj1OAUc1vdr2DD2upvzvIavjMueq2k2Zae0ZmT6RLwte4J7PZvCn0yaH9BI3fvhPvi98jc6Ok3jpgnuaPqBdD/jlm7DxE1gyw7pK+NsnrZEpyZ2sDv3SPKj2Xgmcmgn9L7AuFOsaxLoPSh0mNClEUs36zKGMPgLoPpytMdbQtH4de4cYVBN2fGdd1Tt6GiS2rbfr0bOncsorX/Hu9n8yYuMAzjgyuFrLHz5+lv/lP0cHhjB74uOBN7mIwFFnWbfSfFgz0+o8Ls2zEm9SB2u+qOwR1hW1oXZMK3UY0KQQSeFKCskdWJPSGfAwrFv41mU+iDHWMpnJnXxOfJcQG8eMsf/kotkXc+sXN/NW29c4Oj2zWS/x+7lP8cmep2jDsXw06ZngL+NP7gDHXx3csUqpWjpmLZKqvX0Kwa6nUMfK+I4keTwcHeYlOOtZP8ea6fOU2/2Oke6TnsW9wx7G7Sjil7OnsLOoIKBTV1RXMf71W5i399+04zg+uuRZXQ1NqRZAk0IkubyTh4U6+gjY4HCSXV2NY/fKkM/lk6sSPv4/SO8DQ65stOiF/U7kmqPvptKxgzFvX8LC7esbLf9j/m5OfeVyNld9Qt/EcXx22X9omxiZC3OUUo3TpBBJ4Rp9BOyijOxqF2z7KuRz+bToGWtis7P+EtBwz5tPPJ9bjn0YtxRyzWcTufDNO/g+d1O9Muvycrn83T8zftZYimQ952XezFsX/yXwaQKUUrbTv8ZICtPoowPlpbicRWSYZPhxvjVEM5xK8qzVz448E444PeDDrhxyBgO7HMEdn/+d9eVzuOqzj4hxdSFWUqg2pVQ7dyPiob1jMH8aeSuje/leS0ApFT2aFCIpTB3NC3esR8TQtu0A+HGO9SGekh6GAL3m328N5TzzL80+dFBmTz6+7N/8sHMLT3z3BhsL11DlKSPV2ZnMpKFMzZlky3oISqnw0KQQSWFKCst2bwSg65HjYdNsa8nHYVNCjc6y8wdY+pI1F1B68BPtDcrsWXvlr1Lq0KF9CpEUptFHGwo2A5Bz7NnWrKkr3wo1MourCmZOtaZmHnVHeM6plDqkaFKIpNqaQmh9CrnF2xB3WzokpVrTO+R+D/k/hh7fF3+Hvath7D8OulBNKXV40KQQSbVJIbTRRwXVP5Ek3rngj5torZv79T9Ci23bt/Dlw9aC8kefHdq5lFKHLE0KkVQdep+Cx+Ohkj2kJ3injE7rYs3ns+w1awWrYJTmw7vXQNsecM5DQcemlDr0aVKIJFe5tVJWCOv9/liwB5zl9EjN/nnjiN9Zq3F9EcQHelkBvDQOyvbBhOchIS3o2JRShz5NCpHkqgy5P2FR7joA+nassy5zWqa1EtkPr1jTRAeqJiHs2wiTXtfZQ5VS9iYFETlbRNaLyCYROWg4i4hcISJ5IrLMe7vGzniirro85JFHq/ZaVwnnZDaYCO/UadbC7jOnWh/2TSnJs1ZSy9sAk16D3qeGFJdSqnWwLSmIiBP4F3AOcAwwSUR8XbX0pjFmoPf2nF3xtAiuipCvUdh8YCvG42RQZoN1FGITrWUly/bBjHFQus/3CYyBjfPg6ZOs9X4nvtasq5aVUq2bnTWFocAmY8xmY0wV8AYw3sbXa/nCkBR2lW8n1pPue76groNh0huQvwmeOQUWPuVd5N1AwWZrPqPpo+DVCdbaxL/+HI7UhKCU+pmdVzR3BXbUeZ4LDPNR7kIRORnYAPzOGLPDR5nWoboi5OajYtdO2sQ0slj9EafB5Fkw72747x3WLT4NKous/el9resQBkwKyxTeSqnWJdrTXMwGXjfGVIrItcAM4KDGbRGZAkwB6N69e2QjDCdXeUg1BWsivH10TR7ZeMFuQ+Gq/8KuFdbaxIW5kNEPeo6CjkcE/fpKqdbPzqTwE9CtzvMs77Zaxpj8Ok+fA3yOqTTGTAemA+Tk5JjwhhlBrsqQksKCLSsR8TCwU4ATynU5zroppVSA7OxT+B44UkR6ikgcMBGYVbeAiHSp83QcsNbGeKKvujykabO/3WEtqHNy9oBwRaSUUvXYVlMwxrhE5EbgY8AJvGCMWS0i9wGLjTGzgN+KyDjABRQAV9gVT4vgqghpiot1BeswnnhyumoTkFLKHrb2KRhj5gBzGmy7u87jO4HDZ35lV0VIF6/tKv+RJLoR43SGMSillPqZXtEcSZUlEJcU1KEut5syyaVLYq+mCyulVJA0KUSKxw3lBZAc3Appi3/ahDgq6duhT9OFlVIqSJoUIqUsH4wHkjsFdfgXW5cDcEKWrmuslLKPJoVIKc2z7oNcS3nJnhUY42R0Lx1iqpSyT7QvXjt8lOy17oNsPtpSvIYEutMmIbg+CaWUCoTWFCKlZoK6IJqPSisrKZMt9EgO8KI1pZQKkiaFSCmtqSl0bPahczcsRhwuhnUZHOaglFKqPk0KkVKaB45YSGzX7EM/3/odAOcePTzcUSmlVD2aFCKlJM/qTxBp9qFrClYgrvb0y+jWdGGllAqBJoVIKc0LqunI4/GQ715PRpxen6CUsp8mhUgp3Qspze9knr3ue3AWc3znoTYEpZRS9WlSiJSa5qNmen31hxjjYMqQ82wISiml6tOkEAnGeJuPmpcUPB4Pa4u+Js30Ibt9cFdCK6VUc2hSiITKInBXNrv56JNNy/DE5HFil9E2BaaUUvVpUoiE2gvXmldTeGH5uxgj/HrIOBuCUkqpg2lSiIQgprjYWrCXNSUf08lxPEenZ9oUmFJK1adJIRJqJsNrRlK46/MnQaq468Tf2hSUUkodzNakICJni8h6EdkkInf42B8vIm969y8SkWw744maA9ut+wD7FL7Zto4VxXPo5Die04/Q9ZiVUpFjW1IQESfwL+Ac4Bhgkog0nNHtamC/MeYI4DHgQbviiRpXJSx6BjIHQUpGk8W3H8jjhk+nIsbJX0cdlEeVUspWdk6dPRTYZIzZDCAibwDjgTV1yowH7vU+fgd4UkTEGGNsjCuylvwHCrfDuH/6neLC4/GwZm8ur62cx4fbZ+BxFHHn4H8wrPuRkY1VKXXYszMpdAV21HmeCwzzV8YY4xKRQqADsC/cwbwy90Fey3253jbj4zM6kGzkq4yvbYLBiYfKrB7sW/BnWHCfj+M8eByliKMKgDjpxu8G/YlLB4wKIBKllAqvQ2KRHRGZAkwB6N69e1DnSE3qQGdPysHn9vV6ARQSH1lAGhQ0IlRJAnlx3ejkSPB5MhEhNbYNmSmZnNRtIBP6nUSM0+n336GUUnayMyn8BNSd1jPLu81XmVwRiQHaAPkNT2SMmQ5MB8jJyQmqaWn8Kdcw/pRrgjlUKaUOG3aOPvoeOFJEeopIHDARmNWgzCxgsvfxBODzVtWfoJRShxjbagrePoIbgY8BJ/CCMWa1iNwHLDbGzAKeB14WkU1AAVbiUEopFSW29ikYY+YAcxpsu7vO4wrgIjtjUEopFTi9olkppVQtTQpKKaVqaVJQSilVS5OCUkqpWpoUlFJK1ZJD7bIAEckDtgVxaEdsmD4jDFpqXKCxBaOlxgUaWzBaalzQ/Nh6GGOanL//kEsKwRKRxcaYnGjH0VBLjQs0tmC01LhAYwtGS40L7ItNm4+UUkrV0qSglFKq1uGUFKZHOwA/WmpcoLEFo6XGBRpbMFpqXGBTbIdNn4JSSqmmHU41BaWUUk1odUlBRM4WkfUisklEDlrkWETiReRN7/5FIpIdgZi6ich8EVkjIqtF5CYfZUaJSKGILPPe7vZ1Lpvi2yoiK72vu9jHfhGRx73v2QoRGRyhuI6u834sE5EiEbm5QZmIvG8i8oKI7BWRVXW2tReReSKy0Xvfzs+xk71lNorIZF9lbIjt7yKyzvvzel9E2vo5ttGfvU2x3SsiP9X5mY3xc2yjf8s2xPVmnZi2isgyP8fa/Z75/LyI2O+bMabV3LCm6P4R6AXEAcuBYxqUuQF42vt4IvBmBOLqAgz2Pk4FNviIaxTwYZTet61Ax0b2jwHmYi0bdwKwKEo/291YY60j/r4BJwODgVV1tj0E3OF9fAfwoI/j2gObvfftvI/bRSC2M4EY7+MHfcUWyM/eptjuBW4N4Ofd6N9yuONqsP8R4O4ovWc+Py8i9fvW2moKQ4FNxpjNxpgq4A1gfIMy44EZ3sfvAKeJiK9VOcPGGLPLGLPU+7gYWIu1PvWhYjzwkrEsBNqKSJcIx3Aa8KMxJpgLF0NmjPkCa82Puur+Ls0Azvdx6FnAPGNMgTFmPzAPONvu2IwxnxhjXN6nC7FWPow4P+9bIAL5W7YlLu/nwcXA6+F6veZo5PMiIr9vrS0pdAV21Hmey8EfvrVlvH80hUCHiEQHeJurBgGLfOweLiLLRWSuiPSLVEyAAT4RkSVirYfdUCDvq90m4v+PNFrvW4YxZpf38W4gw0eZlvDeXYVV0/OlqZ+9XW70Nm294KcZJJrv20hgjzFmo5/9EXvPGnxeROT3rbUlhRZNRFKAd4GbjTFFDXYvxWoaGQA8AcyMYGgjjDGDgXOAqSJycgRfu0liLec6Dnjbx+5ovm+1jFV3b3FD+UTk/wAX8KqfItH42T8F9AYGAruwmmpakkk0XkuIyHvW2OeFnb9vrS0p/AR0q/M8y7vNZxkRiQHaAPl2ByYisVg/4FeNMe813G+MKTLGlHgfzwFiRaSj3XF5X+8n7/1e4H2sqntdgbyvdjoHWGqM2dNwRzTfN2BPTTOa936vjzJRe+9E5ApgLHCp90PkIAH87MPOGLPHGOM2xniAZ/28ZlTeN+9nwi+AN/2VicR75ufzIiK/b60tKXwPHCkiPb3fLicCsxqUmQXU9MhPAD739wcTLt42yueBtcaYR/2U6VzTtyEiQ7F+NpFIVskiklrzGKuDclWDYrOAX4nlBKCwTjU2Evx+c4vW++ZV93dpMvCBjzIfA2eKSDtvM8mZ3m22EpGzgduBccaYMj9lAvnZ2xFb3f6oC/y8ZiB/y3Y4HVhnjMn1tTMS71kjnxeR+X2zqwc9WjeskTIbsEYu/J93231YfxwACVjNEJuA74BeEYhpBFZVbwWwzHsbA1wHXOctcyOwGmuUxULgxAi9X728r7nc+/o171nd2AT4l/c9XQnkRPDnmYz1Id+mzraIv29YSWkXUI3VTns1Vl/UZ8BG4FOgvbdsDvBcnWOv8v6+bQKujFBsm7Dalmt+32pG3GUCcxr72Ucgtpe9v0crsD7oujSMzfv8oL9lO+Pybv9Pze9WnbKRfs/8fV5E5PdNr2hWSilVq7U1HymllAqBJgWllFK1NCkopZSqpUlBKaVULU0KSimlamlSUEopVUuTglJKqVqaFJSyiYg8ISJLReT4aMeiVKA0KShlA+8UCJ2Aa7HmH1LqkKBJQakQiUiiiPxPRJw124wxpViLpSwAHveWixORL7yTrinVImlSUKoZvJMCNvy7uQp4zxjjrlOuA5AEFGNNXY2xFov5DLgkQuEq1WyaFJRqgohke9cKfglrRsxuDYpcysEzVk4DHsaaNK3uwj8zveWVapE0KSgVmCOBfxtj+pk6S4J6p3XuZYzZWmdbNnAi1pz8a6mfFFYB2vGsWixNCkoFZpux1qduqCNwoMG2+4H7jDUFcb2k4G1iqqqZk1+plkY7vJQKTKmf7eVYa3QAICIDsVbuGiEi//LuW9ngmHigwo4glQqV1hSUCoExZj/gFJGaxPAg1oJO2caYbGAAdWoK3g7ofcaY6ogHq1QAtKagVOg+waoZeIAkY8ynNTuMMXtEJEVE2htjCoDRwEfRClSppujKa0qFSEQGA78zxlweQNn3gDuMMRvsj0yp5tPmI6VCZIxZCsyve/GaL96RSjM1IaiWTGsKSimlamlNQSmlVC1NCkoppWppUlBKKVVLk4JSSqlamhSUUkrV0qSglFKqliYFpZRStf4fXV0gGLxyUEQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f31bead3510>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "dist = output['RDF']['distance']\n",
    "NaNa = output['RDF']['Na-Na']\n",
    "NaCl = output['RDF']['Na-Cl']\n",
    "ClCl = output['RDF']['Cl-Cl']\n",
    "\n",
    "plt.xlabel(r'r ($\\AA$)')\n",
    "plt.ylabel('g(r)')\n",
    "plt.plot(dist,NaNa,label='Na-Na')\n",
    "plt.plot(dist,NaCl,label='Na-Cl')\n",
    "plt.plot(dist,NaNa,label='Cl-Cl')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation of Ion Pair Lifetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPL distance calculation 14.07% complete"
     ]
    }
   ],
   "source": [
    "ip.runionpair(comx,comy,comz,Lx,Ly,Lz,moltypel,moltype,tsjump,dt,output,ver,IPL_skip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ion Pair Lifetime Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chlorine around sodium: 0.974829130494\n",
      "Sodium around Clorine: 0.967411999742\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFF9JREFUeJzt3X/wXXV95/HnKwkJqCgg0UECDUKcNuuoMN9Fp7pdu2stMDZpZ5cOjDvayso4LTt2dX/guIMuO7tT66y1TGm7dNaxulaWumrjEhe10nGnK0hQwAQmECMK0Uo0LGXLjxDy3j/uyenl8r05yZec3Hy/5/mY+c733M/5fO99f77nfvPKOZ97zklVIUkSwLJZFyBJOnYYCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWqtmHUBh+vUU0+ttWvXzroMSVpUbr/99h9X1equfosuFNauXcuWLVtmXYYkLSpJvnco/Tx8JElqGQqSpJahIElqGQqSpJahIElq9RYKST6W5KEkW6esT5JrkuxIcleS8/qqRZJ0aPrcU/g4cMFB1l8IrGu+Lgf+sMdaJEmHoLdQqKqvAXsO0mUj8IkauQU4KclpfdVz2/17+MiXtrN33/6+XkKSFr1ZzimcDjww9vjBpu1ZklyeZEuSLbt3717Qi33zew9zzVd3sG+/oSBJ0yyKieaquq6q5qpqbvXqzrO0JUkLNMtQ2AWcMfZ4TdMmSZqRWYbCJuBtzaeQXgc8UlU/nGE9kjR4vV0QL8mngTcCpyZ5EPgAcBxAVf0RsBm4CNgBPAb8el+1jKs6Gq8iSYtTb6FQVZd2rC/gN/t6/UnJ0XolSVq8FsVEsyTp6DAUJEktQ0GS1BpcKDjPLEnTDSYUgjPNktRlMKEgSepmKEiSWoaCJKk1uFAoT2mWpKkGEwqe0SxJ3QYTCpKkboaCJKllKEiSWoaCJKk1uFDws0eSNN3gQkGSNJ2hIElqGQqSpJahIElqDS4UvMqFJE03mFCI17mQpE6DCQVJUjdDQZLUMhQkSa3hhYITzZI01WBCwWlmSeo2mFCQJHUzFCRJLUNBktQaXCiUM82SNFWvoZDkgiTbk+xIcuU8689McnOSbyW5K8lF/dXS1zNL0tLRWygkWQ5cC1wIrAcuTbJ+otu/A26oqnOBS4A/6KseSVK3PvcUzgd2VNXOqtoLXA9snOhTwAub5RcBP+ixHklShxU9PvfpwANjjx8EXjvR54PAl5L8C+D5wJt6rEeS1GHWE82XAh+vqjXARcAnkzyrpiSXJ9mSZMvu3buf0wt66WxJmq7PUNgFnDH2eE3TNu4y4AaAqvo6cDxw6uQTVdV1VTVXVXOrV69eUDHOM0tStz5D4TZgXZKzkqxkNJG8aaLP94F/DJDkZxiFwnPbFZAkLVhvoVBV+4ArgJuAexh9ymhbkquTbGi6vRd4Z5I7gU8Dv1blAR5JmpU+J5qpqs3A5om2q8aW7wZe32cNkqRDN+uJ5qPO3RBJmm4woeA9miWp22BCQZLUzVCQJLUMBUlSy1CQJLUGFwqeBiFJ0w0mFPzwkSR1G0woSJK6GQqSpJahIElqDS4UnGaWpOkGEwrOM0tSt8GEgiSpm6EgSWoZCpKk1uBCwROaJWm64YSCpzRLUqfhhIIkqZOhIElqGQqSpNbgQqE8p1mSphpMKDjNLEndBhMKkqRuhoIkqWUoSJJawwsF55klaarBhIInNEtSt8GEgiSpm6EgSWoZCpKkVq+hkOSCJNuT7Ehy5ZQ+v5rk7iTbkvxpn/WA88ySdDAr+nriJMuBa4FfAB4EbkuyqaruHuuzDngf8PqqejjJS3qrx3OaJalTn3sK5wM7qmpnVe0Frgc2TvR5J3BtVT0MUFUP9ViPJKlDn6FwOvDA2OMHm7ZxrwBekeSvktyS5IL5nijJ5Um2JNmye/funsqVJM16onkFsA54I3Ap8MdJTprsVFXXVdVcVc2tXr36KJcoScPRZyjsAs4Ye7ymaRv3ILCpqp6qqu8C9zIKCUnSDPQZCrcB65KclWQlcAmwaaLP5xntJZDkVEaHk3b2WBPlx48kaareQqGq9gFXADcB9wA3VNW2JFcn2dB0uwn4SZK7gZuBf11VP+mjHi9zIUndevtIKkBVbQY2T7RdNbZcwHuaL0nSjM16olmSdAwxFCRJrcGFQnmhC0maajCh4DyzJHUbTChIkroZCpKklqEgSWoNLhQ8o1mSpus8eS3J8cBbgH8AvAx4HNgK3FhV2/ot78jxjGZJ6nbQUEjy7xkFwl8CtwIPAcczukbRbzeB8d6quqvnOiVJR0HXnsI3quoDU9Z9pLlT2plHuCZJ0owcdE6hqm4ESHLx5LokF1fVQ1W1pa/iJElH16FONL/vENuOec4zS9J0XXMKFwIXAacnuWZs1QuBfX0WdqTFc5olqVPXnMIPgNuBDc33Ax4F/mVfRUmSZuOgoVBVdwJ3JvlUVT11lGqSJM3IQecUknwhyS9NWffy5i5q7+inNEnS0dZ1+OidjO6K9tEke4DdjM5TWAt8B/j9qvrzXis8wspTmiVpqq5QeAGwqar+TZK1wGmMzmh+IbCrqr7Tb3lHkPPMktSp6yOpHwUeAaiq+6vq61V1B/Bws06StIR0hcJLq+rbk41N29peKpIkzUxXKJx0kHUnHMlCJEmz1xUKW5K8c7IxyT/nmectLBrOM0vSdF0Tzb8FfC7JW/m7EJgDVgK/0mdhR5rzzJLUrevktR8BP5vk54FXNs03VtVXe69MknTUdd5kB6CqbgZu7rkWSdKMDe52nJKk6QwFSVLLUJAktQYTComfP5KkLr2GQpILkmxPsiPJlQfp90+SVJK5PuuRJB1cb6GQZDlwLXAhsB64NMn6efqdCLwbuLWvWiRJh6bPPYXzgR1VtbOq9gLXAxvn6fcfgA8BT/RYiyTpEPQZCqcDD4w9frBpayU5Dzijqm7ssY5n8DIXkjTdzCaakywDPgK89xD6Xp5kS5Itu3fvXtjrLeinJGlY+gyFXcAZY4/XNG0HnMjo0hl/meR+4HXApvkmm6vquqqaq6q51atX91iyJA1bn6FwG7AuyVlJVgKXAJsOrKyqR6rq1KpaW1VrgVuADVW1pceaJEkH0VsoVNU+4ArgJuAe4Iaq2pbk6iQb+npdSdLCHdIF8RaqqjYDmyfarprS94191tK+Ds40S9I0AzqjedYVSNKxbzChIEnqZihIklqGgiSpNbhQ8IxmSZpuMKHgRLMkdRtMKEiSuhkKkqSWoSBJag0uFJxnlqTpBhMK8eLZktRpMKEgSepmKEiSWoaCJKk1uFAoT2mWpKkGEwqe0SxJ3QYTCpKkboaCJKllKEiSWoMLBaeZJWm6wYWCJGk6Q0GS1DIUJEktQ0GS1DIUJEmtwYWCV7mQpOkGEwrxOheS1GkwoSBJ6mYoSJJahoIkqTXAUHCmWZKm6TUUklyQZHuSHUmunGf9e5LcneSuJH+R5Kd6q6WvJ5akJaS3UEiyHLgWuBBYD1yaZP1Et28Bc1X1KuAzwO/0VY8kqVufewrnAzuqamdV7QWuBzaOd6iqm6vqsebhLcCaHuuRJHXoMxROBx4Ye/xg0zbNZcAXe6xHktRhxawLAEjyz4A54B9OWX85cDnAmWeeucDXGH3f7zyzJE3V557CLuCMscdrmrZnSPIm4P3Ahqp6cr4nqqrrqmququZWr169oGKWNamw3+tcSNJUfYbCbcC6JGclWQlcAmwa75DkXOC/MAqEh3qsheXLRqHwtLsKkjRVb6FQVfuAK4CbgHuAG6pqW5Krk2xoun0YeAHwZ0nuSLJpytM9Z8tjKEhSl17nFKpqM7B5ou2qseU39fn649xTkKRugzmjedky5xQkqctgQuHvDh/NuBBJOoYNJhSWNSP18JEkTTeYUFjRpIKHjyRpusGEwvJmpPvcU5CkqQYTCu3Ja4aCJE01mFDwI6mS1G0woXBgT+Fp5xQkaarBhMKK5R4+kqQugwmFA+cpONEsSdMNJhQ8o1mSug0mFLwgniR1G04o+OkjSeo0mFDw8JEkdRtMKKxY5kSzJHUZTCh4RrMkdRtMKBzXnKfw1NOGgiRNM5hQWLViOQBP7vOGCpI0zYBCYTTUJ556esaVSNKxazChsGxZWLliGU/sMxQkaZrBhALA8SuW8eRTHj6SpGmGFQrHLedJ9xQkaapBhcKq45bxhHsKkjTVoELh+BXLnWiWpIMYVCicsHI5jxsKkjTVoELhRSccxyOPPzXrMiTpmDWoUDjl+SvZ87d7Z12GJB2zBhUKJz9vJXv+n6EgSdMMKhRWn7iKR5/cx2N79826FEk6Jg0qFM55yQsA2P7Xj864Ekk6Ng0qFNaf9kIAtv3gb2ZciSQdm3oNhSQXJNmeZEeSK+dZvyrJf2/W35pkbZ/1rDn5BM445QQ+961d7Hvak9gkaVJvoZBkOXAtcCGwHrg0yfqJbpcBD1fVOcDvAh/qq56mJn7jjedw+/ce5l3/7Zts3fUIj+/1vAVJOmBFj899PrCjqnYCJLke2AjcPdZnI/DBZvkzwO8nSVV/N1K+9PwzeXzv0/zHzffwlXt+BMDzVi7neStXsGrFMpLRXdqWZRQi6Xi+yULn61+T6zJanneQz+o8vY+3C1qcut5Tz/0HDoNvokXl3W9ax8bXnN7ra/QZCqcDD4w9fhB47bQ+VbUvySPAi4Efj3dKcjlwOcCZZ575nAt7xxvO4i2vPo1bd+7h+3seY8/f7uWxvfvYu68oiip4en894++lqkjm/+s80Dr59zX5Mwce7x/LvMP5e6+x/geet89/L3TkTfs3ePz/QZPvmfHt3odp72sde178/FW9v0afoXDEVNV1wHUAc3NzR+T/Ni858Xh+6dUvOxJPJUlLRp8TzbuAM8Yer2na5u2TZAXwIuAnPdYkSTqIPkPhNmBdkrOSrAQuATZN9NkEvL1Z/qfAV/ucT5AkHVxvh4+aOYIrgJuA5cDHqmpbkquBLVW1CfivwCeT7AD2MAoOSdKM9DqnUFWbgc0TbVeNLT8BXNxnDZKkQzeoM5olSQdnKEiSWoaCJKllKEiSWllsnwBNshv43gJ//FQmzpYeAMc8DI55GJ7LmH+qqlZ3dVp0ofBcJNlSVXOzruNocszD4JiH4WiM2cNHkqSWoSBJag0tFK6bdQEz4JiHwTEPQ+9jHtScgiTp4Ia2pyBJOojBhELX/aIXsyT3J/l2kjuSbGnaTkny5ST3Nd9PbtqT5Jrm93BXkvNmW/2hSfKxJA8l2TrWdthjTPL2pv99Sd4+32sdC6aM94NJdjXb+Y4kF42te18z3u1JfnGsfdG875OckeTmJHcn2Zbk3U37Ut7O08Y8u21dVUv+i9FVWr8DvBxYCdwJrJ91XUdwfPcDp060/Q5wZbN8JfChZvki4IuMbub1OuDWWdd/iGP8OeA8YOtCxwicAuxsvp/cLJ8867Edxng/CPyrefqub97Tq4Czmvf68sX2vgdOA85rlk8E7m3GtpS387Qxz2xbD2VPob1fdFXtBQ7cL3op2wj8SbP8J8Avj7V/okZuAU5KctosCjwcVfU1RpdXH3e4Y/xF4MtVtaeqHga+DFzQf/WHb8p4p9kIXF9VT1bVd4EdjN7zi+p9X1U/rKpvNsuPAvcwumXvUt7O08Y8Te/beiihMN/9ovu9+/XRVcCXktze3M8a4KVV9cNm+a+BlzbLS+l3cbhjXApjv6I5VPKxA4dRWILjTbIWOBe4lYFs54kxw4y29VBCYal7Q1WdB1wI/GaSnxtfWaP9ziX9MbMhjBH4Q+Bs4DXAD4H/PNty+pHkBcD/AH6rqv5mfN1S3c7zjHlm23oooXAo94tetKpqV/P9IeBzjHYlf3TgsFDz/aGm+1L6XRzuGBf12KvqR1X1dFXtB/6Y0XaGJTTeJMcx+sfxU1X12aZ5SW/n+cY8y209lFA4lPtFL0pJnp/kxAPLwJuBrTzz/tdvB/68Wd4EvK355MbrgEfGds0Xm8Md403Am5Oc3OyOv7lpWxQm5n5+hdF2htF4L0myKslZwDrgGyyy932SMLpF7z1V9ZGxVUt2O08b80y39axn34/WF6NPKtzLaIb+/bOu5wiO6+WMPmlwJ7DtwNiAFwN/AdwHfAU4pWkPcG3ze/g2MDfrMRziOD/NaDf6KUbHSy9byBiBdzCanNsB/Pqsx3WY4/1kM567mj/408b6v78Z73bgwrH2RfO+B97A6NDQXcAdzddFS3w7TxvzzLa1ZzRLklpDOXwkSToEhoIkqWUoSJJahoIkqWUoSJJahoIGJ8lJSX5j7PHLknymp9f65SRXLeDnVib5WpIVfdQlTeNHUjU4zTVm/mdVvfIovNb/ATZU1Y8X8LMfYHSRs08d+cqk+bmnoCH6beDs5jr1H06yNs19C5L8WpLPN9ftvz/JFUnek+RbSW5JckrT7+wk/6u5COH/TvLTky+S5BXAkwcCIcnHk/xRki1J7k3ylqb97yX5RlPPXUnWNU/xeeCtR+MXIh3grqmG6ErglVX1Gmj3HMa9ktHVKo9ndEbsv62qc5P8LvA24KOM7pX7rqq6L8lrgT8A/tHE87we+OZE21pG17E5G7g5yTnAu4Dfq6pPNZcoWN703Qr8/ec2VOnwGArSs91co2vbP5rkEeALTfu3gVc1V7T8WeDPRpeuAUY3PZl0GrB7ou2GGl3k7L4kO4GfBr4OvD/JGuCzVXUfQFU9nWRvkhObeqTeefhIerYnx5b3jz3ez+g/UsuA/1tVrxn7+pl5nudxRnsb4yYn8aqq/hTY0PTfnGR8j2MV8MQCxyEdNkNBQ/Qoo1sfLkiNrnf/3SQXQ3uv4FfP0/Ue4JyJtouTLEtyNqOLGW5P8nJgZ1Vdw+gKoK9qnvfFwI+r6qmF1iodLkNBg1NVPwH+KsnWJB9e4NO8FbgsyYGr085368OvAedm7BgT8H1Glzr+IqM5iSeAXwW2JrmD0XzGJ5q+Pw/cuMD6pAXxI6lSj5L8HvCFqvpKko8z+ijsIZ0TkeSzjG5Yf2+fNUrj3FOQ+vWfgOcd7g81n0L6vIGgo809BUlSyz0FSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktf4/nA5AsGlLKeAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f31b99a8c50>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "Correlation = output['Ion_Pair_Lifetime']['Na around Cl correlation']\n",
    "Time = output['Ion_Pair_Lifetime']['Correlation_Time']\n",
    "plt.xlabel('time (ps)')\n",
    "plt.ylabel('C(t)')\n",
    "plt.plot(Time,Correlation)\n",
    "\n",
    "print('Chlorine around sodium: {}'.format(output['Ion_Pair_Lifetime']['Cl around Na'][4]))\n",
    "print('Sodium around Clorine: {}'.format(output['Ion_Pair_Lifetime']['Na around Cl'][4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation of MSD, Diffusivity and Nernst-Einstein Conductivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "unwrap complete\n",
      "MSD calculation 49.97% complete"
     ]
    }
   ],
   "source": [
    "output = m.runMSD(comx, comy, comz, Lx, Ly, Lz, Lx2, Ly2, Lz2, moltype, moltypel, dt, tsjump, output, ver, MSD_skip, MSD_num_init)\n",
    "cd.calcdiffusivity(output, moltypel, dt,D_Tolerance)\n",
    "output = ne.calcNEconductivity(output, molcharge, Lx, Ly, Lz, nummoltype, moltypel, temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diffusivity Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diffusivity of sodium: 1.32568724574e-08 m^2/s\n",
      "Diffusivity of chlorine: 1.20292297281e-08 m^2/s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEKCAYAAAAmfuNnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd4VVXWwOHfAkLvRVoIoSO9RIoFRYqABbvgCFg+0VFn1HGk2GiioKIjigUVuyhDZ6QLWFDEoJBCDT2UUAIJLaSt749z0CuSkITbkqz3ee6Te/dpa3tjFuecffYSVcUYY4zxtSKBDsAYY0zhYAnHGGOMX1jCMcYY4xeWcIwxxviFJRxjjDF+YQnHGGOMX1jCMcYY4xeWcIwxxviFJRxjjDF+USzQAQRa1apVNTw8PNBhGGNMvrJmzZpDqlotN9sU+oQTHh5OZGRkoMMwxph8RUR25nYbu6RmjDHGLyzhGGOM8Qu/JRwRqSMiy0VkvYjEisijbntlEVkiIlvcn5XcdhGRiSISJyJRItLOY1+D3PW3iMggj/b2IhLtbjNRRMRf/TPGGJM9f97DSQeeUNVfRaQcsEZElgB3A9+o6jgRGQYMA4YCvYFG7qsj8DbQUUQqAyOACEDd/cxV1SPuOvcDPwPzgV7AgtwGmpaWRnx8PCkpKRfU4UApWbIkoaGhhISEBDoUY4z5nd8SjqruA/a574+JyAagNtAXuMpd7WNgBU7C6Qt8ok7BnlUiUlFEarrrLlHVRAA3afUSkRVAeVVd5bZ/AtxIHhJOfHw85cqVIzw8nPx2kqSqHD58mPj4eOrVqxfocIwx5ncBuYcjIuFAW5wzkepuMgLYD1R339cGdntsFu+2Zdcef472XEtJSaFKlSr5LtkAiAhVqlTJt2dnxpiCy+8JR0TKAjOAx1Q12XOZezbj8xKkIjJYRCJFJPLgwYNZrePrMHwmP8dujCm4/JpwRCQEJ9l8rqoz3eYE91IZ7s8DbvseoI7H5qFuW3btoedo/wtVnayqEaoaUa1arp5bMsaYfG/P0VOMnree9IxMvx7Xn6PUBPgA2KCqr3osmgucGWk2CJjj0T7QHa3WCUhyL70tAnqKSCV3RFtPYJG7LFlEOrnHGuixr3xHRHjiiSd+//zKK68wcuTIwAVkjMn3MjOVT1ftpOer3/LlL7tYvy/5/Bt5kT/PcC4DBgBXi8ha99UHGAf0EJEtQHf3MzijzLYBccB7wEMA7mCBMcAv7mv0mQEE7jrvu9tsJQ8DBoJFiRIlmDlzJocOHQp0KMaYAmDHoRP0f28Vz86OoW1YJRY91oVWoRX9GoM/R6n9AGR1c6HbOdZX4OEs9jUFmHKO9kigxQWEGTSKFSvG4MGDee211xg7duyfls2bN4/nn3+e1NRUqlSpwueff0716tWz2JMxpjDLyFSm/LCdCUs2EVK0CONvacntEXUCcq+30M+ldj6j5sWyfq93Tzub1SrPiOubn3e9hx9+mFatWjFkyJA/tV9++eWsWrUKEeH999/npZdeYsKECV6N0RiT/21OOMaQ6VGs3X2U7hdfxNibWlK9fMmAxWMJJ4iVL1+egQMHMnHiREqVKvV7e3x8PHfccQf79u0jNTXVnrcxxvxJWkYm76zYyhvL4ihbshiv92vDDa1rBXwEqyWc88jJmYgvPfbYY7Rr14577rnn97Z//OMf/Otf/+KGG25gxYoVNpjAGPO7mD1JPDk9ig37krm+dS1GXt+MKmVLBDoswCbvDHqVK1fm9ttv54MPPvi9LSkpidq1nWdaP/7440CFZowJIilpGby0cCN9J63k8PHTTB7Qnjf6tw2aZAOWcPKFJ5544k+j1UaOHMltt91G+/btqVq1agAjM8YEgzU7E7l24ve8tWIrt7SrzZLHr6Rn8xqBDusv7JJakDp+/Pjv76tXr87Jkyd//9y3b1/69u0biLCMMUHkZGo6Ly/axEc/7qBWhVJ8cm8HujQO3ofZLeEYY0w+tDLuEMNmRrE78RQDO9dlSK+mlC0R3H/Sgzs6Y4wxf5KcksaL8zcwdfVu6lUtw1eDO9GxfpVAh5UjlnCMMSaf+GZDAk/PiuHAsRQe6FKfx3s0pmRI0UCHlWOWcIwxJsgdOZHKqHmxzF67lybVy/HugPa0ruPfaWm8wRKOMcYEKVVlfvR+RsyN4ejJNB7t1oiHuzakeLH8OcDYEo4xxgShA8dSeG52LAtj99OydgU+va8jF9csH+iwLkj+TJOFwP79++nXrx8NGjSgffv29OnTh82bN9OiRYGYm9QYkwVVZdovu+k+4VuWbTrAsN5NmfXQpfk+2YCd4QQlVeWmm25i0KBBfPnllwCsW7eOhISEAEdmjPGlnYdPMHxmND9uPUyHepV58eaWNKhW1jcHO7wVqjTwzb6zYAknCC1fvpyQkBAefPDB39tat27Njh07AheUMcZn0jMy+eCH7by2dDMhRYow9qYW9L8kjCJFfDDZ5tHdsPgZ2DAXHvgOarT0/jGyYAnnfBYMg/3R3t1njZbQe1yWi2NiYmjfvr13j2mMCUqxe5MYNiOa6D1J9GhWnTF9W1Cjgg9KCKSlwE9vwPevgipc9RRUaeT942TDbwlHRKYA1wEHVLWF2/YV0MRdpSJwVFXbiEg4sAHY5C5bpaoPutu0Bz4CSuFUBX1UVVVEKgNfAeHADuB2VT3i844ZY0wepKRl8Po3W5j83TYqlQ5h0p3t6NOyhm9KCGxaCAuHwpEdcPENcM1YqBjm/eOchz/PcD4C3gQ+OdOgqneceS8iE4Akj/W3qmqbc+znbeB+4GechNMLp5T0MOAbVR0nIsPcz0MvOOpszkR8pXnz5kyfPt3vxzXG+MeqbYcZPjOa7YdOcFv7UJ6+9mIqli7u/QMd3goLh8GWxVC1CQyYDQ26ev84OeS3UWqq+h2QeK5l4qT024Gp2e1DRGoC5VV1lVuC+hPgRndxX+DMXP0fe7TnO1dffTWnT59m8uTJv7dFRUWxe/fuAEZljLlQySlpDJ8ZTb/Jq0jPzOSz+zry8m2tvZ9sUk/A0lHwVifY+RP0HAt/XxnQZAPBcw/nCiBBVbd4tNUTkd+AZOAZVf0eqA3Ee6wT77YBVFfVfe77/UB1H8fsMyLCrFmzeOyxxxg/fjwlS5YkPDyc//znP4EOzRiTR4tj9/PsnBgOHjvN/VfU4/EejSld3Mt/glUhdiYsfhaS90Dr/tB9FJQLjj+HwZJw+vPns5t9QJiqHnbv2cwWkRyX3nTv6WhWy0VkMDAYICzM/9cxc6JWrVpMmzbtL+0xMTEBiMYYk1cHjqUwcm4s86P307RGOSYPiPDNtDQJsbBgKOz4Hmq0gls/hLCO3j/OBQh4whGRYsDNwO/DslT1NHDafb9GRLYCjYE9QKjH5qFuG0CCiNRU1X3upbcDWR1TVScDkwEiIiKyTEzGGJNXqsp/I+N5/uv1pKRn8uQ1TRjcpT4hRb18J+PUUVjxIqx+D0qWh+teg3aDoEjwTeoZ8IQDdAc2qurvl8pEpBqQqKoZIlIfaARsU9VEEUkWkU44gwYGAm+4m80FBgHj3J9z/NkJY4w5Y+fhEzw1K5qVcYfpEF6ZF2/xwQOcmZmw9nNYOhJOJUL7e+DqZ6B0Ze8ex4v8OSx6KnAVUFVE4oERqvoB0I+/DhboAowWkTQgE3hQVc8MOHiIP4ZFL3Bf4CSaaSJyH7ATZxBCnqmqb4Yn+oEznsIY42/pGZlMWbmdV5f4+AHOPWtg/pPOzzodoc9MqNnau8fwAb8lHFXtn0X73edomwHMyGL9SOAvE4qp6mGg24VF6ShZsiSHDx+mSpUq+S7pqCqHDx+mZEkfPDhmjMnS+r3JDJsZRVR8Et0vvogxN7agZoVS3j3IiUPOGc1vn0HZi+Cmd6HVHZBP/k4FwyW1oBMaGkp8fDwHDx4MdCh5UrJkSUJDQ8+/ojHmgqWkZfDGsi28++02KpYO4c0723Jty5re/cdqRjpEfgDLxzpDni99BLoMce7Z5COWcM4hJCSEevXqBToMY0yQ+2VHIkNnRLHt4AlublebZ69tRqUyXn6mZscPMH8IHIiF+l2h93io1uT82wUhSzjGGJNLx1LSeGnhJj5dtZPaFUvx8b0duLJxNe8eJGkPLHkWYmZAhTC44zNoel2+uXx2LpZwjDEmF5ZtTODpWTHsT07hnsvC+XfPJpQp4cU/pemn4adJ8N0rkJkOVw6Fyx6D4qW9d4wAsYRjjDE5cPj4aUbNW8/cdXtpXL0sk/52Ke3CKnn3IFuWOA9vJm51zmauGQuVwr17jACyhGOMMdlQVWav3cPoees5fjqdx7o34qGrGlK8mBcf4EzcBgufgs0LoEpDuGsGNOzuvf0HCUs4xhiThfgjJ3lmdgwrNh2kbVhFxt/SisbVy3nvAKkn4IfXYOVEKFLMmfes00NQzAczRwcBSzjGGHOWzEzlk5928NIipyTXiOubMbBzOEW99QDn2ZNstrwdeoyC8rW8s/8gZQnHGGM8bEk4xtAZUfy66yhdGldj7I0tqFPZizfs90c792l2rnQm2bzlA6jb2Xv7D2KWcIwxBkhNz+Sdb7fy5rI4Spcoyqu3t+amtrW99wDnyUTnwc3IKVCyIlz3H2g3MCgn2fQVSzjGmEJv7e6jDJ0exaaEY1zXqiYjb2hO1bIlvLPzzAxY8xEsGwMpSXDJ/8FVw4N6kk1fsYRjjCm0TqamM2HxZj5cuZ2LypXk/YERdG/mxWJlO390ZglIiIbwK6DXOKjxl6kgCw1LOMaYQumHLYcYPiuK3YmnuKtTGEN7NaVcyRDv7Dx5rzMgIGY6lA+F2z6CZjfm61kCvMESjjGmUDl6MpXnv97A9DXx1K9ahmkPdKZDPS9d3kpLgVWT4LsJziwBXYbA5Y9B8TLe2X8+ZwnHGFMoqCrzo/czYm4sR06m8nDXBvzj6kaUDPHCTXtV2LwQFg6HI9sL5CwB3mAJxxhT4CUkp/DM7BiWrE+gRe3yfHzvJTSvVcE7Oz+0BRYOg7ilULUJDJgFDa72zr4LGC8X186aiEwRkQMiEuPRNlJE9ojIWvfVx2PZcBGJE5FNInKNR3svty1ORIZ5tNcTkZ/d9q9EpGA+qmuMyTFVZerqXXSf8C3fbT7I8N5Nmf3QZd5JNinJzn2atzrD7tVwzQvw95WWbLLhzzOcj4A3gU/Oan9NVV/xbBCRZjilp5sDtYClItLYXTwJ6AHEA7+IyFxVXQ+Md/f1pYi8A9wHvO2rzhhjgtuuwycZOiOKn7YdpnP9Krx4c0vCq3rhXkpmJkR9BUtHwPEEaHsXdBvhVOA02fJnienvRCQ8h6v3Bb5U1dPAdhGJAzq4y+JUdRuAiHwJ9BWRDcDVwJ3uOh8DI7GEY0yhk5GpfPzjDl5etIliRYQXb25Jv0vqeOcBzj2/OrMExK+G2u2h31QIbX/h+y0kguEeziMiMhCIBJ5Q1SNAbWCVxzrxbhvA7rPaOwJVgKOqmn6O9Y0xhUTcgeMMmb6OX3cd5eqmFzH2phbUrFDqwnd8/CAsGw2/fgplqkLft6B1fyjit7sSBUKgE87bwBhA3Z8TgHt9fVARGQwMBggLC/P14YwxPpaekcm7323j9W+2ULp4Uf5zRxv6tql14Wc1GWnwy/uw/EVIOwGdH4Yrh0BJLw04KGQCmnBUNeHMexF5D/if+3EPUMdj1VC3jSzaDwMVRaSYe5bjuf65jjsZmAwQERGhF9gNY0wArd+bzJAZ64jZk0yfljUYdUMLqpXzwrQ021bAgmFwcIMzEKDXeKjW+LybmawFNOGISE1V3ed+vAk4M4JtLvCFiLyKM2igEbAaEKCRiNTDSSj9gDtVVUVkOXAr8CUwCJjjv54YY/ztdHoGk5bF8daKrVQsHcLbf2tH75Y1L3zHR3fD4qdh/RyoWBf6fQFN+hT6WQK8wW8JR0SmAlcBVUUkHhgBXCUibXAuqe0AHgBQ1VgRmQasB9KBh1U1w93PI8AioCgwRVVj3UMMBb4UkeeB34AP/NQ1Y4yfrd19lCHT17E54Tg3t63Ns9c1o1KZC3wSIi0FfnwDvp/gfO76NFz6TwgpeeEBGwBEtXBfUYqIiNDIyMhAh2GMyYGUtAxeW7KZ977fRvXyJXnhppZ0beqF4cibFsLCoXBkB1x8gzNLQEW7v5sdEVmjqhG52SbQgwaMMSZHVm9PZOiMKLYfOkH/DmEM79OU8hc62ebhrc50NFsWQdXGMGA2NOjqnYDNX1jCMcYEtROn03lp4UY+/mkndSqX4vP/68hlDate2E5TTziXzn58A4qWgJ7PQ4cHoJhNUOJLlnCMMUHr+y0HGTYjmr1Jp7jnsnCevKYJpYtfwJ8tVVg/GxY9Dcl7oFU/6DEKytXwXtAmS5ZwjDFBJ+lUGmO/Xs+0yHjqVyvDfx/oTET4BZYQOLAB5j8JO76H6i3hlg+gbmfvBGxyxBKOMSaoLFmfwDOzozl0PJW/X9WAR7tdYAmBlCRYMR5+fgdKlIM+r0DEvVDEC2UJTK5YwjHGBIXEE6mMmhfLnLV7aVqjHO8PvISWoRfwRP+ZSTaXPAcnDkL7QXD1c1CmiveCNrliCccYE1CqytfR+xgxJ5bklDQe796Yv1/VgOLFLmCesr1rnctn8auhdgTc+RXUbue9oE2eWMIxxgTMgeQUnp0Tw6LYBFqFVuDzWzvStEb5vO/wZCIsGwORH9okm0HIEo4xxu9UlRm/7mH0vFhS0jMZ1rsp/3d5PYoVzWNiyMyANR85ySYlGTo+CFcNg1IVvRq3uTCWcIwxfrXn6CmemhnNt5sPElG3EuNvbUWDamXzvsNdP8P8f8P+KKh7OfR5Gao3817Axmss4Rhj/CIzU/li9S5enL8BBUbd0JwBnepSpEgeJ8U8lgBLR8K6L6BcLbh1CjS/2SbZDGKWcIwxPrfz8AmGzohi1bZELmtYhXE3t6JO5dJ521lGGqyeDCvGQdopuPxxuOLfUOICzpKMX1jCMcb4TEam8uHK7byyeBMhRYow/paW3B5xAeWet30LC4bAwY3QsLtTo6ZqQ+8GbXzGEo4xxifiDhzjyelR/LbrKN2aXsTzF1LuOSnemY5m/Wy3Rs1UaNLbLp/lM5ZwjDFelZaRyeTvtvH60i2ULnGB5Z7TT/9Ro0Yz3Ro1/4CQPCYuE1CWcIwxXhO7N4kh06OI3ZvMtS1rMvKG5nkv97x5sVOjJnEbXHw99BwLlep6N2DjV/6s+DkFuA44oKot3LaXgeuBVGArcI+qHhWRcGADsMndfJWqPuhu0x74CCgFzAcedUtMVwa+AsJxqoferqpH/NE3Ywq70+kZvLksjrdXbKVi6eK8c1c7erXIY7nnxG2w8CnYvACqNIK7ZkLDbt4N2ASEPx+//QjodVbbEqCFqrYCNgPDPZZtVdU27utBj/a3gfuBRu7rzD6HAd+oaiPgG/ezMcbHftt1hOsm/sAby+K4oU0tlv6rS96STepJWPY8TOrkzOjcYwz8/UdLNgWI385wVPU798zFs22xx8dVwK3Z7UNEagLlVXWV+/kT4EZgAdAXuMpd9WNgBTD0wiM3xpzLqdQMXl2yiQ9+2E718iX58J5L6NokD+WeVWH9HFj8DCTthpa3Q4/RUD6PZ0gmaAXTPZx7cS6JnVFPRH4DkoFnVPV7oDYQ77FOvNsGUF1V97nv9wPVfRyvMYXWz9sOM3RGFDsOn+TOjmEM792Ucnkp93xgozPMefu3UL0F3DwZ6l7q/YBNUAiKhCMiTwPpwOdu0z4gTFUPu/dsZotI85zuz72no9kcbzAwGCAsLCzvgRtTyBw/nc74BRv5dNVOwiqX5ov7O3JpgzyUe05Jhm/dGjXFy0Dvl50aNUWD4k+S8ZGAf7sicjfOYIJuqqoAqnoaOO2+XyMiW4HGwB4g1GPzULcNIEFEaqrqPvfS24Gsjqmqk4HJABEREVkmJmPMH77bfJDhM51yz/deVo9/X9M49+WeVSFqGix5Fo4fgHYDoNsIZ2ZnU+AFNOGISC9gCHClqp70aK8GJKpqhojUxxkcsE1VE0UkWUQ6AT8DA4E33M3mAoOAce7POX7sijEFlme55wbVyjD9wc60r5uHcs/7opwaNbtXQe320H+q89MUGv4cFj0V56Z+VRGJB0bgjEorASxxHwo7M/y5CzBaRNKATOBBVU10d/UQfwyLXuC+wEk000TkPmAncLsfumVMgfbNhgSemnWB5Z5PJsLysRA5BUpVhhvehDZ/sxo1hZC4V7EKrYiICI2MjAx0GMYElSNuuefZbrnnl29tnftyz5kZ8Osn8M1oSDkKl9wPXYdDqUq+Cdr4lYisUdWI3GwT8Hs4xpjgsjBmH8/MjuXoyVQe7daIh7s2zH25592/ODVq9q2FupdB75egRgvfBGzyDUs4xhgADh0/zYg5sXwdvY/mtcrzyb0daFYrl+Wejx+ApaNg7WdQribc8gG0uMUm2TSAJRxjCj1VZV7UPkbOjeV4SjpPXtOEwV3qE5Kbcs8Z6fDLe7D8BadGzWWPQpcnoUQ53wVu8h1LOMYUYgeSU3h6dgxL1ifQuk5FXr61FY2r5zJJbP/eeXjzwHpocLVz+axqI98EbPK18yYcEemBM+JrkqquFZHB7nMsxph8SlWZ+eseRv9vPSlpGTzVpyn3XlaPYrk5q0na40xHEzsTKoTBHZ9D02vt8pnJUk7OcO4F/g48487I3Ma3IRljfGlf0imemhnN8k0HiahbiZdubUX9arkoz5x+Gn6aBN+97NSouXIYXP6Y1agx55WThHNMVY8C/xaRccAlPo7JGOMDqspXv+xm7NcbSM9URlzfjIGdwylaJBdnJFuWOpfPErdC0+vgmrFQKdxnMZuCJScJ5+szb1R1mIj8w4fxGGN8YHfiSYbPjOaHuEN0ql+Z8be0om6VMjnfQeJ2WPQUbJoPVRrCXTOgYXffBWwKpPMmHFWdc9bnN7Ja1xgTXDIzlc9/3sm4BRsBeP7GFtzZIYwiOT2rST0JK/8DP/wHihSD7qOg00NQrLgPozYFlVdGqYlIRfeymzEmSOw8fIIh06P4eXsiVzSqyos3tyS0UumcbawKG+Y5ZzVJu6HFrdBzDJSv5dugTYGWk1Fq7XFmc56IU0Kg+VmvFkAZoKLvwjTG5FRGpvLRjzt4edFGQooW4aVbWnFbRCiS09FjBzc792m2LYeLmsPdX0P45b4N2hQKOTnDeRd4ANgFHANigY3ABqAf0EZVsywFYIzxn60HjzNkehRrdh7h6qYX8cJNLalRoWTONj59DL59CVa9BSFlnOdpIu6zGjXGa3Lym/Qj8CTwK1AaeE9VpwGIyJOWbIwJvPSMTN7/YTuvLtlMqZCivHp7a25qWztnZzWqEP1fWPwsHN8Pbe+CbiOhbDWfx20Kl5wMGviniJRW1ZPuczjPiMjjwGigcE81bUwQ2JxwjCf/u4518Ulc07w6Y25swUXlcnhWsz8a5g+BXT9CrbbQ73MIzdUEwMbkWI7Olc8UR3Nr0vxLROoCzwPVRaSrqi73YYzGmHNIy8jknRVbmbhsC+VKhvBG/7Zc16pmzs5qTh2BZWMh8gOnXMD1E6HtAKtRY3wqTxdnVXUnMEBEJgDjRGSkql7p3dCMMVmJ3ZvEkOlRxO5N5rpWNRl1Q3OqlC1x/g0zM+G3T+GbUU7SueT/oOtTVqPG+MUF/XNGVdeqai9gZE7WF5EpInJARGI82iqLyBIR2eL+rOS2i4hMFJE4EYkSkXYe2wxy198iIoM82tuLSLS7zUTJ8bAcY/KH1PRMXl28ib5vriQh+TTv3NWeN+9sl7NkE78G3u8G8/4JVRvDA99Bn5ct2Ri/8cr5cy4uqX0E9DqrbRjwjao2Ar5xPwP0Bhq5r8HA2+AkKJzy1B2BDsCIM0nKXed+j+3OPpYx+VZU/FGuf+MHJi6L44bWtVj6ry70alHj/BsePwhzHoH3r4bkvXDze3DPAqjR0vdBG+PBr+MdVfU7EQk/q7kvcJX7/mNgBTDUbf9EnRrYq0SkoojUdNdd4t5PQkSWAL1EZAVQXlVXue2fADcCC3zXI2N8LyUtg/8s3cLk77ZSrVwJPhgUQbeLq59/w4x05x7NsrGQdgIu/SdcOcRq1JiAyVPCEZFqAKp60AsxVFfVfe77/cCZ/5NqA7s91ot327Jrjz9HuzH51pqdRxgyfR1bD57gjog6PHXtxVQoFXL+DXeshPlPwoFYqN/VeaamWmPfB2xMNnKccNz7ISOAR3AuxYmIpANvqOpobwSjqioiPh9qLSKDcS7TERYW5uvDGZNrp1IzmLB4Ex+s3E6tCqX45N4OdGmcg+dikvc6z9PETIcKdeCOz5xZne12pgkCuTnDeRy4DLhEVbcDiEh94G0ReVxVX8tjDAkiUlNV97mXzM48SLoHqOOxXqjbtoc/LsGdaV/htoeeY/2/cAvITQaIiIiwZ4lMUPl522GGzohix+GT3NUpjGG9L6ZsifP8r5qe6swQ8O1LkJkOVw6Fyx6D4jmcO80YP8jNoIEBQP8zyQZAVbcBdwEDLyCGucCZkWaDgDke7QPd0WqdgCT30tsioKeIVHIHC/QEFrnLkkWkk3s2NtBjX8YEvROn0xkxJ4Y7Jq8iU+GL+zvy/I0tz59s4pbC251h6QiofyU8/LMz1NmSjQkyuTnDCVHVQ2c3qupBEcnBRWUQkak4ZydVRSQe5xLdOGCaiNwH7MQpZw0wH+gDxAEngXvc4yWKyBjgF3e90WcGEAAP4YyEK4UzWMAGDJh84cethxg6I4r4I6e4+9JwhvRqQuni5/nf88gOWPQ0bPwfVG4Af5sOjXr4JV5j8iI3CSc1j8t+p6r9s1jU7RzrKvBwFvuZAkw5R3skzuzVxuQLx0+nM27BBj5btYvwKqX5anBnOtSrnP1Gaadg5evww2sgRaDbCOj8MBTLwbM4xgRQbhJOaxFJ9vjseRcyhxM3GWPOWBl3iCHTo9ibdIr7Lq/Hv3s2oVTxollvoAobv4ZFw+HoLmh+M/R8HirYYEyTP+Q44ahqNv8nGGNy6lhKGi8u2MgXP++iftUyTH+wM+3rnues5tAWWDAUtn4D1S6GQfNBgyDwAAAdG0lEQVSgXhf/BGyMl+RmWPQlwG5V3e9+HgjcAuwARnncRzHGZOH7LQcZNiOafUmnGNylPv/q0ZiSIdn8W+70cfjuZfhpEoSUgl7jnPnPiubotqkxQSU3l9TeBboDiEgXnJv9/wDa4AwxvtXr0RlTQCSnpPHi/A1MXb2bBtXKMP3vl9IuLJs5zFQhZgYsfgaO7YM2d0H3EVD2Iv8FbYyX5SbhFPU4i7kDmKyqM4AZIrLW+6EZUzCs2HSA4TOjSUhO4cErG/BY90bZn9Xsj3FKPO9cCTXbwO2fQp1L/BewMT6Sq4QjIsVUNR1nVNngPO7HmEIh6VQaY79ez7TIeBpeVJaZD11GmzoVs97g1FFY/gL88h6UrADX/QfaDYQidvvUFAy5SRRTgW9F5BBwCvgeQEQaAkk+iM2YfGv5Rues5sCxFB66qgH/7JbNWU1mJqz9HJaOhFOJ0P4euPoZKH2egQTG5DO5GaU2VkS+AWoAi93nZMAZHv2IL4IzJr9JOpnG6P+tZ8av8TSuXpZ3B1xG6+zOavascSbZ3LMG6nSEPjOhZmv/BWyMH+VmlNpcj4/3etQ2E0CBG7wYlzH5ztL1CTw1K5rDJ1J5pGtD/tGtISWKZXFWc+KQU3Xz10+dgQA3vQut7rBJNk2BlptLap1xygJMBX7mzw9+GlNoHT2Zyuh565n52x6a1ijHB4MuoWVohXOvnJEOaz6EZWMg9YQzQ8CVQ6Fkef8GbUwA5Cbh1AB6AP2BO4GvgamqGuuLwIzJD5a4ZzVHTqTyz26NeKRrQ4oXy2JO3J0/OpfPEmKg3pVOjZqLmvo3YGMCKDf3cDKAhcBCESmBk3hWiMgoVX3TVwEaE4yOnEhl1LxYZq/dy8U1y/Ph3ZfQonYWZzXJ+2DJcxA9DcqHwu2fwMU32OUzU+jkajizm2iuxUk24cBEYJb3wzImeC2M2c8zs2M4ejKVx7o34qGrsjirSU+Fn9+Bb8dDRip0eRIufxyKl/F/0MYEgdwMGvgEZybm+ThT2cT4LCpjglDiiVRGzI1l3rq9NKtZnk/u7UCzWlnce9m6zJn77NBmaNwLrnkBqjTwb8DGBJncnOHcBZwAHgX+efYoNVW1u56mwJofvY9nZ8eQnJLGEz0a8+BVDQgpeo6zmqO7YNFTsGEeVKoHd06Dxtf4P2BjglBu7uHkpjqoMQXCoeOnGTEnlq+j99Gidnk+v60jTWuc499Waadg5UT44VWnRs3Vz0LnRyDEKncYc0bAp6QRkSbAVx5N9YHngIrA/cBBt/0pVZ3vbjMcuA/IAP6pqovc9l7A60BR4H1VHeeXTpgC6X9Re3luTizHU9J58pomDO5S/69nNaqwaQEsHAZHd0Lzm9waNaGBCdqYIBbwhKOqm3BmnEZEigJ7cAYi3AO8pqqveK4vIs2AfkBzoBawVEQau4sn4Qzdjgd+EZG5qrreLx0xBcaBYyk8NzuWhbH7aR1agZdva03j6uX+uuLhrc59mrglUK0pDJwL9a/0f8DG5BMBTzhn6QZsVdWdkvWQ0b7Al6p6GtguInFAB3dZnKpuAxCRL911LeGYHFFVZvy6hzH/W8+ptAyG9mrK/VfUo9jZZzWnj8P3rzg1aoqWcAYEdBhsNWqMOY9gSzj9cGYyOOMRt9BbJPCEqh4BagOrPNaJd9vAmQnBs72jD2M1Bcieo6d4amY0324+SETdSoy/tRUNqpX980qqEDsTFj0Dx/ZC6zuh+0goVz0QIRuT7wRNwhGR4jjzsQ13m94GxuDM0zYGmADc66VjDcYtrxAWFuaNXZp8KjNT+Xz1LsbN34ACI69vxsDO4RQpctYZdsJ6p0bNju+hRiu47SMIs3/PGJMbQZNwgN7Ar6qaAHDmJ4CIvAf8z/24B6jjsV2o20Y27X+iqpNxqpQSERGh51rHFHw7Dp1gyIwoVm9P5PKGVXnx5pbUqVz6zyudOgorxsHqyc58Z9e9Bu0GWY0aY/IgmBJOfzwup4lITVXd5368CTjzoOlc4AsReRVn0EAjYDXO80CNRKQeTqLphzPnmzF/kpGpTPlhOxOWbCKkaBHG39KS2yPq8Kf7hpmZsG4qLB3hzOwccY8z1Nlq1BiTZ0GRcESkDM7osgc8ml8SkTY4l9R2nFmmqrEiMg1nMEA68LA7zxsi8giwCGdY9BSbWNScbXPCMZ6cHsW63UfpfvFFjL2pJdXLn/WszN7fnEk243+B0A7wt+lQq01gAjamAJE/6qgVThERERoZGRnoMIyPpaZn8vaKrby5fAvlSoYw4vpm3NC61p/Pak4chmWjYc3HUKYa9Bjt1KgpYs88G3M2EVmjqhG52SYoznCM8aXo+CSenL6OjfuPcX3rWoy8vhlVypb4Y4XMDIicAsueh9PHoNNDcNVQKJnF7M/GmDyxhGMKrJS0DF7/ZguTv9tGlTLFeW9gBD2anTWEedcqmP9v2B8N4VdAn5fhoosDE7AxBZwlHFMgRe5IZMiMKLYdPMHtEaE8fW0zKpTyeDDz2H5YMgKivoTyteHWD51paaxGjTE+YwnHFCgnTqfz8qJNfPzTDmpVKMUn93agS+Nqf6yQkQY/v+sMdc44DVc84bysRo0xPmcJxxQYP2w5xLCZUcQfOcWgznUZ0qspZUp4/IpvXe7WqNkEjXpCr3FWo8YYP7KEY/K9pFNpvPD1Br6K3E29qmWY9kBnOtTzeF7m6G5Y/DSsnwOVwqH/V9CkV8DiNaawsoRj8rWl6xN4enY0B4+d5oEr6/N498aUDHFnAUhLgZ/egO8mOJ+7PgOX/sNq1BgTIJZwTL50+PhpRs1bz9x1e2laoxzvDYygVWjFP1bYtBAWDoUjO6BZX+g5FirWyXJ/xhjfs4Rj8hVVZV7UPkbOjeVYShqPdW/EQ1c1pHgx9+HMw1th4XDYsgiqNoYBs6FB18AGbYwBLOGYfCQhOYVnZsewZH0CrUMr8NKtnWhSwy2MlnoCvp8AP77h1Kjp+Tx0eACKFQ9s0MaY31nCMUFPVflvZDxjvl5Panomw3s35b7L3cJoqrB+Nix6GpL3QKt+0GMUlKsR6LCNMWexhGOC2u7Ekzw1K5rvtxyiQ3hlxt3SkvpnCqMd2AgLnoTt30GNlnDrFAjrFNiAjTFZsoRjglJmpvLpqp2MX7gRAcb0bc7fOtZ1CqOlJMGK8bD6XSheFq6dAO3vsRo1xgQ5Szgm6Gw7eJyhM6L4ZccRrmjkFEYLrVTaqVGz9ktY8hycOAjtB8HVz0GZKoEO2RiTA5ZwTNBIz8jk/R+289qSzZQoVoSXbm3Fbe1DnRICe9e6NWpWQ+0IuPMrqN0u0CEbY3LBEo4JChv3JzNkehRR8Un0bFad529swUXlS8LJRFg2BiI/hDJVoe9b0Lq/1agxJh8KmoQjIjuAY0AGkK6qESJSGfgKCMep+nm7qh4Rp2rW60Af4CRwt6r+6u5nEPCMu9vnVfVjf/bD5E5qeiaTlsfx1oo4ypcM4c0723Jty5qIZjo1ar4ZDSnJ0PFBuGoYlKp4/p0aY4JS0CQcV1dVPeTxeRjwjaqOE5Fh7uehQG+gkfvqCLwNdHQT1AggAqc09RoRmauqR/zZCZMzv+06wrAZ0WxKOEbfNrUYcX1zKpcpDrtXOzVq9q2Dupc7NWqqNwt0uMaYCxRsCedsfYGr3PcfAytwEk5f4BN16mOvEpGKIlLTXXeJqiYCiMgSoBcw1b9hm+ycOJ3OK4s38dGPO6heriTvD4yge7PqcCwBZo2EdV9AuVrOMOfmN1uNGmMKiGBKOAosFhEF3lXVyUB1Vd3nLt8PnCnXWBvY7bFtvNuWVbsJEis2HeDpWTHsOXqKAZ3qMqRXE8qFAD9NcmrUpJ2Cyx+HK/4NJcoGOlxjjBcFU8K5XFX3iMhFwBIR2ei5UFXVTUYXTEQGA4MBwsLCvLFLcx6Hj59mzP/WM3vtXhpUK8N/H+zMJeGVYdu3sGAIHNwIDbtDr/FQtWGgwzXG+EDQJBxV3eP+PCAis4AOQIKI1FTVfe4lswPu6nsAz6l/Q922PfxxCe5M+4pzHGsyMBkgIiLCK0nMnJuqMnvtHkbPW8/x0+n8s1sjHu7agBIn9sF//wWxs6BiXeg3FZr0tstnxhRgQZFwRKQMUERVj7nvewKjgbnAIGCc+3OOu8lc4BER+RJn0ECSm5QWAS+ISCV3vZ7AcD92xXjYnXiSp2fH8N3mg7QNq8i4m1vRpGpx+PE1Z6JNzYSuT7s1akoFOlxjjI8FRcLBuTczyxntTDHgC1VdKCK/ANNE5D5gJ3C7u/58nCHRcTjDou8BUNVEERkD/OKuN/rMAALjPxmZyocrtzNh8WaKCIy6oTl3dapL0bglMG0oJG6Di693atRUqhvocI0xfiLOQK/CKyIiQiMjIwMdRoGxYV8yw2ZEsS4+ia5NqvH8TS2pnbkPFj4FmxdAlUbQezw07BboUI0xF0BE1qhqRG62CZYzHJPPpaRl8MayLbz77TYqlAphYv+2XH9xBeSHCbByIhQNgR5jnAc4rUaNMYWSJRxzwVZtO8xTM6PZdugEt7QL5Zk+Tam0cwFMegaSdkPL26HHaChfM9ChGmMCyBKOybOkU2mMW7CBqat3U6dyKT69rwNXVEyEGbfB9m+hegu4eTLUvTTQoRpjgoAlHJMnC2P289ycGA4dP839V9Tj8S41KP3jBPj5HSheBnq/DBH3QlH7FTPGOOyvgcmVhOQURsyJZWHsfi6uWZ73B7anVeJiePdGOH4A2g2AbiOcmZ2NMcaDJRyTI5mZyleRu3lh/gZS0zMZ2qsp/9foOCEL+8HuVVC7PfSf6vw0xphzsIRjzivuwHGemhXN6u2JdKpfmfF96lB33Wvw/hQoVRlueBPa/M1q1BhjsmUJx2QpNT2Td77dypvL4igZUoSXbm7GbUVWIJ/fCSlHocNguGq41agxxuSIJRxzTmt2JjJsRjRbDhznulY1Gd3+JJVX3A371kLdy6D3S1CjRaDDNMbkI5ZwzJ8kp6Tx8sJNfPbzTmqWL8mn/cK5YsdbMPUzKFcTbvkAWtxik2waY3LNEo753aJYZ6jzgWOnuadzKEMr/0CJBXc7NWouexS6PAklygU6TGNMPmUJx5CQnMJzc2JYFJtA0xrl+Lx7Gg0jH4Rf10ODq53LZ1UbBTpMY0w+ZwmnEMvMVL5YvYvxCzaSmpHJqKsqMuDY+xSZPxMqhsEdn0PTa+3ymTHGKyzhFFJbEo4xfGY0kTuP0KV+OV6v+yOVIieCZjgjzy571GrUGGO8yhJOIXM6PYO3lm/lrRVxlClRjE+6JHNF3LPIT1uh6XVwzVioFB7oMI0xBZAlnEJk9fZEhs+MYuvBE9zbDIbKB5RYvRCqNIS7ZkDD7oEO0RhTgAX80XARqSMiy0VkvYjEisijbvtIEdkjImvdVx+PbYaLSJyIbBKRazzae7ltcSIyLBD9CUZJp9IYPjOa29/9CU09xfKIH3lu5z2U2PU9dB8Ff//Jko0xxueC4QwnHXhCVX8VkXLAGhFZ4i57TVVf8VxZRJoB/YDmQC1gqYg0dhdPAnoA8cAvIjJXVdf7pRdBSFVZGLOfEXNjOXQ8hZeb7+KWQ29RJGY3tLgVeo6B8rUCHaYxppAIeMJR1X3APvf9MRHZANTOZpO+wJeqehrYLiJxQAd3WZyqbgMQkS/ddQtlwtmXdIpnZ8eydEMC11RP4pWLvqDc1u/houZw99cQfnmgQzTGFDIBTzieRCQcaAv8DFwGPCIiA4FInLOgIzjJaJXHZvH8kaB2n9Xe0cchB52MTOWzVTt5edEmSmSeYHaj5bTe8wWSUsZ5nibiPqtRY4wJiKD5yyMiZYEZwGOqmiwibwNjAHV/TgDu9dKxBgODAcLCwryxy6Cwaf8xhs2M4rddRxhaK4r7T39Esd0J0PYu6DYSylYLdIjGmEIsKBKOiITgJJvPVXUmgKomeCx/D/if+3EPUMdj81C3jWza/0RVJwOTASIiItQLXQiolLQM3lwWxzvfbqV9iXgia31J1cQ1UKst9P8CQiMCHaIxxgQ+4YiIAB8AG1T1VY/2mu79HYCbgBj3/VzgCxF5FWfQQCNgNSBAIxGph5No+gF3+qcXgbNq22GemhnNoUMJfFhzAZcfnYOkVILrJ0LbAVajxhgTNAKecHDu1QwAokVkrdv2FNBfRNrgXFLbATwAoKqxIjINZzBAOvCwqmYAiMgjwCKgKDBFVWP92RF/SjqZxgvzNzAtcicPlv+Rf1WYSsjRJLjk/6DrU1CqUqBDNMaYPxHVfH9F6YJERERoZGRkoMPIMVXl6+h9jJy7njqn1jOp4lRqnVgPYZ2hz8tQo2WgQzTGFAIiskZVc3W9PhjOcEwO7Tl6iudmx7B24xZerDCLniGLQWrAze9By9tskk1jTFCzhJMPZGQqn/y0g9cWred2FvNW2ekUT0uBS/8JVw6xGjXGmHzBEk6Q27AvmWEzoykR/xNfl/2MOmnbIayr80xNtcbn34ExxgQJSzhBKiUtg4nfbGHOd5E8U2IqvUv8gJauA70+c2Z1tstnxph8xhJOEPox7hDPzfyVbkkzWVZyNsUlEy4filz2GBQvHejwjDEmTyzhBJEjJ1J5Yf4GEn77miklPiUsZC806gPXvACV6wU6PGOMuSCWcIKAqjJ33V7en7ucf6R/SM/ikWRWagC9p0OjHoEOzxhjvMISToDtTjzJ6FlraLbtQ6aHzCOkeFG4cgRFOj8MxUoEOjxjjPEaSzgBkp6RyUcrt/Pbki8YUeRjQkMOktn8For0HAMVsqvOYIwx+ZMlnACI3ZvEG9MW0P/wJP6vaBRpVZrAdR9SpN4VgQ7NGGN8xhKOH51KzWDS4rWU/fk/vFF0PpQoiXZ7kZAO90PRkECHZ4wxPmUJx0++33yAZdPf5oHTH1Kj6BFSW95J8WtGQdmLAh2aMcb4hSUcH0s8kcoH0+fRZetLjCiykeNVWsBN/6V4nUsCHZoxxviVJRwfUVW+/mUDyQvG8HjmQtKKlyOtx6uUveRuKFI00OEZY4zfWcLxgd2Hj7Po81e58fB7VJLjJLcYQKVrR0LpyoEOzRhjAsYSjhelZ2Qyb8HX1P9lJP8ncRyo1AZue51KtdsEOjRjjAm4ApdwRKQX8DpO1c/3VXWcP467IW4bu/47lL4pS0guVokj3d/gok4DbJJNY4xxFaiEIyJFgUlADyAe+EVE5qrqel8d81TKab6b+hKddrxNI0lhR+O7qXfLKKRkBV8d0hhj8qUClXCADkCcqm4DEJEvgb6ATxLOupULKL10GNfoDraUbU+NOyZSP6yFLw5ljDH5XkFLOLWB3R6f44GOvjjQ6okD6JA4lwSpyuYrJ9H4qr/Z5TNjjMlGQUs4OSIig4HBAGFhYXnaR2bFcH4qdQ9t7xxN9TLlvRmeMcYUSAUt4ewB6nh8DnXb/kRVJwOTASIiIjQvB+o0cExeNjPGmEKrSKAD8LJfgEYiUk9EigP9gLkBjskYYwwF7AxHVdNF5BFgEc6w6CmqGhvgsIwxxlDAEg6Aqs4H5gc6DmOMMX9W0C6pGWOMCVKWcIwxxviFJRxjjDF+YQnHGGOMX1jCMcYY4xeimqfnHgsMETkI7Mzj5lWBQ14MJz+wPhcO1ueC70L7W1dVq+Vmg0KfcC6EiESqakSg4/An63PhYH0u+ALRX7ukZowxxi8s4RhjjPELSzgXZnKgAwgA63PhYH0u+PzeX7uHY4wxxi/sDMcYY4xfWMLJIxHpJSKbRCRORIYFOp6cEJEdIhItImtFJNJtqywiS0Rki/uzktsuIjLR7V+UiLTz2M8gd/0tIjLIo729u/84d1vJ7hg+6uMUETkgIjEebQHrY3bH8HGfR4rIHve7XisifTyWDXfj2SQi13i0n/N32i338bPb/pVb+gMRKeF+jnOXh5/vGF7qbx0RWS4i60UkVkQeddsL7PecTZ/z1/esqvbK5Qun9MFWoD5QHFgHNAt0XDmIewdQ9ay2l4Bh7vthwHj3fR9gASBAJ+Bnt70ysM39Wcl9X8ldttpdV9xte2d3DB/1sQvQDogJhj5mdQw/9Hkk8O9zrNvM/X0tAdRzf4+LZvc7DUwD+rnv3wH+7r5/CHjHfd8P+Cq7Y3ixvzWBdu77csBm95gF9nvOps/56nsO+B/B/PgCOgOLPD4PB4YHOq4cxL2DvyacTUBN931NYJP7/l2g/9nrAf2Bdz3a33XbagIbPdp/Xy+rY/iwn+H8+Y9vwPqY1TH80Oes/hD96XcVp3ZU56x+p3H+gB4Cip39u39mW/d9MXc9yeoYPvy+5wA9CsP3fI4+56vv2S6p5U1tYLfH53i3LdgpsFhE1ojIYLetuqruc9/vB6q777PqY3bt8edoz+4Y/hLIPgbyd+UR9/LOFPnjMmZu+1wFOKqq6We1/2lf7vIkd32/9dm9vNMW+JlC8j2f1WfIR9+zJZzC5XJVbQf0Bh4WkS6eC9X5Z4pPhy364xiBPn6g++h6G2gAtAH2ARMCG473iUhZYAbwmKomey4rqN/zOfqcr75nSzh5sweo4/E51G0Laqq6x/15AJgFdAASRKQmgPvzgLt6Vn3Mrj30HO1kcwx/CWQfA/K7oqoJqpqhqpnAezjfdXbxZNV+GKgoIsXOav/TvtzlFdz1fd5nEQnB+cP7uarOdJsL9Pd8rj7nt+/ZEk7e/AI0ckd1FMe5kTY3wDFlS0TKiEi5M++BnkAMTtxnRucMwrk2jNs+0B190wlIci8lLAJ6ikgl9/S9J8613n1Asoh0ckf0DDxrX+c6hr8Eso9ZHcOnzvxRdN2E812fiaefO/KoHtAI5wb5OX+n3X/FLwdudbc/u29n+nwrsMxdP6tjeKtvAnwAbFDVVz0WFdjvOas+57vv2Vc3tQr6C2dUymackRlPBzqeHMRbH2dEyTog9kzMONdivwG2AEuBym67AJPc/kUDER77uheIc1/3eLRHuL/wW4E3+ePB4nMew0f9nIpzaSEN55ryfYHsY3bH8HGfP3WPF+X+Yajpsf7TbjybcEdfZfc77f7urHb/W/wXKOG2l3Q/x7nL65/vGF7q7+U4l7KigLXuq09B/p6z6XO++p5tpgFjjDF+YZfUjDHG+IUlHGOMMX5hCccYY4xfWMIxxhjjF5ZwjDHG+IUlHGOMMX5hCccYLxKRiiLykPu+lohM9+K+HxORge77puJMR/+biDTIYv0vRaSRt45vzIWy53CM8SJ3YsX/qWoLL++3GPArzhT16W4dk2Kq+nw221wJ3KWq93szFmPyqtj5VzHG5MI4oIGIrMV5Gv1iVW0hIncDNwJlcKYAeQWnHskA4DTQR1UT3bOVSUA14CRwv6puBK4GfnWTTR/gMSBDRLoB1+HUMgnFqXcyRlW/Ar4HPhKRYvrHLMDGBIxdUjPGu4YBW1W1DfDkWctaADcDlwBjgZOq2hb4CWe+LoDJwD9UtT3wb+Att/0yYA2Aqs7HKZD1mqp2BXoBe1W1tXtmtdBdLxNnOpLWvuioMbllCccY/1muqsdU9SBOTZF5bns0EO5OPX8p8F/3DOldnCJfuD8PZrHfaKCHiIwXkStUNclj2QGglrc7Ykxe2CU1Y/zntMf7TI/PmTj/LxbBKYLV5hzbnsKZRPEvVHWziLTDmZTxeRH5RlVHu4tLutsaE3B2hmOMdx3DqTmfa+oU1NouIreBMyW9iJy5HLYBaHiu7USkFs7luc+Al4F2Hosb88eU9cYElJ3hGONFqnpYRFaKSAxOksitvwFvi8gzQAjwJU5JiQU4U9GfS0vgZRHJxClR8HcAEakOnFLV/XmIwxivs2HRxuQTIjILGKKqW3K4/uNAsqp+4NvIjMkZu6RmTP4xjD8GEeTEUeBjH8ViTK7ZGY4xxhi/sDMcY4wxfmEJxxhjjF9YwjHGGOMXlnCMMcb4hSUcY4wxfvH/LVn2zXIMnz4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f31b997d090>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "Na_MSD = output['MSD']['Na']\n",
    "Cl_MSD = output['MSD']['Cl']\n",
    "MSD_time = output['MSD']['time']\n",
    "plt.xlabel('time(fs)')\n",
    "plt.ylabel(r'MSD $\\AA^2$')\n",
    "plt.plot(MSD_time,Na_MSD,label='Na')\n",
    "plt.plot(MSD_time,Cl_MSD,label='Cl')\n",
    "plt.legend()\n",
    "\n",
    "print('Diffusivity of sodium: {} m^2/s'.format(output['Diffusivity']['Na']))\n",
    "print('Diffusivity of chlorine: {} m^2/s'.format(output['Diffusivity']['Cl']))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
